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relativistic  terms,  the  leading  zonal  harmonics  of  earth  and  moon,  the  precession  of 
the  lunar  equator,  and  the  tidal  couple  between  earth  and  moon. 'The  tidal  term  in 
the  moon's  mean  longitude  is  found  to  be  -19"  t  4"  per  century  squared.  The  solution 
also  yields  an  extrapolation  of  the  atomic  time  scale  back  to  1912.5.  At  that  time, 
the  difference  between  atomic  and  ephemeris  time  is  about  6  t  2  seconds.  It  is  found 
that  solar  oblateness  cannot  quite  be  determined  with  optical  data  covering  about  50 
years,  but  J2  is  unlikely  to  be  much  larger  than  lO”^. 

-This  solution  is  believed  to  be  the  only  simultaneous  improvement  of  the  orbits 
of  moon  and  planets.  The  simultaneity  is  found  to  be  an  essential  feature  in 
separating  the  moon's  mean  motion,  the  lunar  tidal  deceleration,  and  the  corrections 
to  the  earth  rotation  rate.  It  is  now  possible  to  refer  all  astronomical  events  of 
the  past  60  years  to  a  time  with  uniform  rate,  namely  the  atomic  clock  system. 
Considering  the  long  baseline,  this  model  should  facilitate  the  prediction  of  fast 
variables,  such  as  the  lunar  longitude,  with  considerably  increased  confidence.  ,  The 
planetary  orbital  elements  compete  with  efforts  of  similar  scope  and  accuracy  at  the 
Massachusetts  Institute  of  Technology  and  the  Jet  Propulsion  Laboratory. 


DD  ,Fr„1473 

'N  0101.807-6801 


(PAGE  |) 


T- 


UNCLASSIFIED 


Security  Classification 


CONTENTS 


Page 

FOREWORD  .  i 

ABSTRACT  .  ii 

I.  INTRODUCTION  .  1 

II.  OBJECTIVES,  RESULTS  OBTAINED  AND  MISSED .  4 

III.  SPECIAL  FEATURES  OF  DESIGN  AND  ANALYSIS  .  7 

IV.  OBSERVATIONS  .  9 

V.  DATA  PREPARATION  .  12 

VI.  ALGORITHM  AND  DERIVATIONS .  16 

1.  PRECOMPUTATIONS .  16 

2.  NUMERICAL  INTEGRATION  ROUTINE  .  17 

3.  EQUATIONS  OF  MOTION . 21 

Spherical  Harmonics  of  Earth  and  Moon  . 22 

Lunar  and  Solar  Tide  Coupling  Effects  . 34 

The  Oblate  Sun  . 37 

4.  INTERPOLATION  AND  LIGHT  TIME  . 39 

5.  RESIDUALS  AND  CORRECTIONS  . 40 

6.  PARTIAL  DERIVATIVES . 42 

Partiais  for  Planetary  Orbits  . 42 

Partials  for  the  Lunar  Orbit  . 45 

Partiais  for  Lunar  Tide  Coupling . 51 

Partiais  for  Solar  Tide  Coupling  . 55 

Partiais  for  Solar  Oblateness  . 57 

Partiais  for  Clock  Corrections . 64 

Partiais  for  Lunar  Declination  Bias . 67 

7.  NORMAL  EQUATIONS  AND  THEIR  SOLUTIONS  . 68 

VII.  CHECKOUT  AND  EXPERIMENTS . 78 

1.  INTEGRATION  ROUTINE  . 78 

2.  TRUNCATION  ERRORS  . 78 

3.  COMPARISON  WITH  OTHER  INTEGRATORS . 82 

4.  PHASE  EFFECTS  FOR  VENUS,  MERCURY,  AND  MOON  ....  84 

5.  OBLATENESS  PERTURBATIONS  . 84 

6.  LUNAR  TESSERAL  HARMONICS . 87 

7.  PARTIALS  FOR  TIDAL  EFFECTS  . 88 

VIII.  CONSTANTS  AND  NUMERICAL  RESULTS  . 89 

IX.  SUMMARY  AND  RECOMMENDATIONS  . 99 

REFERENCES . 100 

Preceding  page  blank 

iii 


CONTENTS  (Continued) 


APPENDICES 

A.  AUXILIARY  TABLES 

VII.  USNO  Transit  Circle  Data  Punched  at  Dahlgren 

VIII.  Observations  on  Dahlgren  Master  Tape 

IX.  Identification  of  Pluto  Observations 

X.  Adams  and  Cowell  Coefficients 

XI.  Number  and  Quality  of  Observations 

XII.  Heliocentric  Equatorial  Coordinates 

XIII.  Lunar  Osculating  Geocentric  Equatorial  Elements 

XIV.  Osculating  Ecliptic  Elements 

XV.  Correlation  Coefficients 

XVI.  Correlation  Coefficients 

XVII.  Correlation  Coefficients 
XVIII.  Correlation  Coefficients 

XIX.  Atomic  Time  1912.5  to  1954.5 

B.  AUXILIARY  FIGURES 

9.  Moon.  Residuals  in  Right  Ascension  Times  cos  5 

10.  Moon.  Residuals  in  Declination 

11.  Sun.  Residuals  in  Right  Ascension  Times  cos  8 

12.  Sun.  Residuals  in  Declination 

13.  Mercury.  Residuals  in  Right  Ascension  Times  cos  8 

14.  Mercury.  Residuals  in  Declination 

15.  Venus.  Residuals  in  Right  Ascension  Times  cos  5 

16.  Venus.  Residuals  in  Declination 

17.  Mars.  Residuals  in  Right  Ascension  Times  cos  5 

18.  Mars.  Residuals  in  Declination 

19.  Jupiter.  Residuals  in  Right  Ascension  Times  cos  8 

20.  Jupiter.  Residuals  in  Declination 

21.  Saturn.  Residuals  in  Right  Ascension  Times  cos  8 

22.  Saturn.  Residuals  in  Declination 

23.  Uranus.  Residuals  in  Right  Ascension  Times  cos  5 

24.  Uranus.  Residuals  in  Declination 

25.  Neptune.  Residuals  in  Right  Ascension  Times  cos  5 

26.  Neptune.  Residuals  in  Declination 

27.  Pluto.  Residuals  in  Right  Ascensioi.  Times  cos  5 

28.  Pluto.  Residuals  in  Declination 

C.  DISTRIBUTION 


t 

! 


iv 


LIST  OF  TABLES 


Table  Page 

I.  Correcting  Catalog  Times  to  UT  . 12 

II.  Maximum  Truncation  Error  of  Moon  After  57  Years .  80 

III.  NWL  vs.  JPL  Coordinates  After  4400  Days .  83 

IV.  Reciprocal  Masses .  90 

V.  Tidal  Coefficient  and  Declination  Bias  .  96 

VI.  Suppressed  Parameters .  98 


v 


LIST  OF  FIGURES 


Figure  Page 

1.  Scheme  for  Generating  Starting  Table  . 19 

2.  Relation  of  Lunar  Equator  to  Ecliptic .  27 

3.  Relation  of  Selenocentric  Longitude  to  Ecliptic .  30 

4.  Relation  of  Argument  of  Latitude  u  to  Polar  and 

Rectangular  Coordinates .  56 

5.  Effective  Truncation  Error  of  Moon  over  55  Years .  81 

6.  Venus  Residuals  in  Right  Ascension  Uncorrected 

for  Phase  Effects .  85 

7.  Effects  on  Lunar  Residuals  due  to  J2  of  Earth  and 

Precession  of  Lunar  Orbit . 7 .  86 

8  Atomic  Minus  Ephemeris  Time,  Extrapolated  Back  to  1912.5  .  .  . .  95 


VI 


I. 


INTRODUCTION 


The  principal  source  of  information  on  the  positions  of  planets  and  moon 
are  the  various  national  ephemerides.  These  publications  enjoy  a  reputation  of  great 
scientific  competence.  Rightfully  so.  Most  of  the  tables  that  are  published  today  are 
based  on  the  life  work  of  many  great  men  in  dynamical  astronomy.  There  is 

abundant  evidence  of  the  painstaking  and  loving  care  that  went  into  various 
theories  of  sun,  moon,  and  planets. 

This  recognition  must  not  impede  our  desire  to  make  progress  when  there  is 
a  need  to  do  so.  And,  in  spite  of  the  renown  of  the  classical  theories,  there  is 
considerable  room  for  improvement.  A  major  reason  for  present  defects  in  the 

ephemerides  lies  in  the  fact  that  a  substantial  part  of  the  theories  was  done  by 

hand,  and  that  quite  some  time  ago.  There  is  only  so  much  algebra  a  man  can  do 
in  forty  years,  so  that  all  series  have  to  be  truncated  somewhere.  Complications 

arise  since  it  is  very  difficult  to  find  all  terms  greater  than  a  given  threshold. 
Moreover,  the  determination  of  constants  in  such  a  theory  cannot  go  much  beyond 
the  orbital  elements  of  the  body  under  investigation.  Although  it  is  no  secret  that 
changes  in  one  constant  require  the  simultaneous  adjustment  of  many  others,  it  is 
simply  not  possible  to  treat  50  or  more  parameters  by  hand.  Such  a  partial  solution 
usually  leads  to  formal  variances  that  are  too  small.  At  times  they  are  sti  bad  that 
completely  unrealistic  .  error  estimates  result.  Finally,  not  every  theory  is  fitted 
directly  to  observations.  No  one  doubts  that  this  is  definitely  the  best  way  to 

proceed,  but  it  is  also  known  to  be  considerably  more  laborious.  Hence,  some 

theories  made  use  of  a  previous  investigator’s  residuals  instead  of  actual  observations. 
It  is  clear  that  this  can  and  will  produce  further  complications. 

As  a  consequence  of  the  deficiencies  just  discussed,  is  there  any  direct 
evidence  for  the  shortcomings  of  the  presently  published  ephemerides?  Yes,  there  are 
some  large  discrepancies  between  observations  and  published  places.  But  it  is  only 
fair  to  note  that  there  are  cases  in  which  agreement  is  still  excellent,  it  is  difficult, 
for  example,  to  find  fault  with  Newcomb’s  motion  of  the  sun,  even  though  the 
theory  is  75  years  old.  We  had  occasion  to  compare  his  position  predictions  for  the 

present  time  with  our  own  solution  and  found  the  differences  to  be  much  less  than 

1".  However,  on  the  other  end  of  the  spectrum  are  Pluto  and  the  moon.  The 
discrepancies  in  Pluto’s  published  places  are  now  about  10".  This  cannot  come  as  a 
surprise  since  the  tables  are  based  on  a  set  of  orbital  elements  determined  in  1932! 
The  residuals  in  the  transit  circle  observations  of  the  moon,  as  published  in  several 
series,  are  also  quhe  large  and  presently  exceed  20".  However,  to  quote  this  number 
may  possibly  be  unfair  to  Brown’s  lunar  theory.  Since  empirical  terms  and  other 
corrections  are  involved,  it  is  undoubtedly  possible  to  explain  or  remove  such  large 
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planning,  the  moral  support  and  encouragement  of  Gerald  M.  Clemence  and  Paul 

Herget  is  gratefully  acknowledged.  At  the  Naval  Observatory,  quite  a  number  of 

people  were  most  helpful.  In  the  Six-Inch  Division,  we  would  like  to  express  special 
thanks  to  David  K.  Scott,  who,  with  his  profound  knowledge  of  practical 

astronomy,  has  set  us  straight  on  seveial  occasions.  B.  L.  Klock  and  A.  N.  Adams 
were  equally  willing  to  listen  to  our  problems  and  to  provide  immediate  assistance. 

Klock  also  was  kind  enough  to  let  us  use  the  newly  reduced  lunar  observations 

before  they  were  released  for  publication.  In  the  Nautical  Almanac  Office, 

Raynor  L.  Duncombe  and  his  staff  assisted  us  on  many  occasions.  They  were 
particularly  helpful  in  locating  much  of  the  valuable  observational  material. 
Thomas  C.  Van  Flandern  brought  his  knowledge  of  lunar  theory  to  bear  on  several 
problems.  Since  Douglas  A.  O’Handley  at  JPL  is  engaged  in  work  similar  to  ours,  it 
was  natural  that  we  had  a  large  number  of  discussions  with  him  on  many  aspect  of 
this  study.  Although  work  of  J.  Derral  Mulholland  on  the  lunar  orbit  follows  a 
somewhat  different  philosophy,  there  was  some  overlap  and  exchange  of  information. 
On  a  few  occasions  we  compared  notes  with  Irwin  I.  Shapiro  of  MIT.  His  keen 
observations  always  provide  much  food  for  thought.  It  should  also  be  remembered 

that  Ash,  Shapiro,  and  Smith  (1967)  were  first  to  publish  a  modern  solution  for 

planetary  orbits  and  masses.  If  any  of  our  friends  were  left  out,  it  was 
unintentional,  and  we  apologize  for  the  oversight. 


II. 


OBJECTIVES,  RESULTS  OBTAINED  AND  MISSED 


Our  original  plans  called  for  a  solution  of  the  orbital  elements  and  masses  of 
the  major  planets  based  on  optical  observations  back  to  1750  or  so  and  on  radar 
ranges  to  Mercury,  Venus,  and  Mars. 

We  soon  found  that  it  was  relatively  easy  to  obtain  meridian  circle 
observations  for  about  the  last  55  years.  Older  data  are  more  difficult  to  resurrect. 
We  decided  to  begin  our  work  with  the  55  year  span  and  to  extend  to  older 
observations  at  a  later  time. 

In  the  early  phases  of  the  study  we  were  not  particularly  interested  in  the 
motion  of  the  moon.  It  was  clear  that  it  had  to  be  considered  but  only  inasmuch 
as  it  affects  the  motion  of  the  earth.  Available  knowledge  of  the  orbit  of  the  moon 
was  sufficient  to  fill  this  need. 

The  next  two  decisions  changed  the  course  of  this  study  drastically.  Since  we 
had  written  an  n-body  program,  it  was  only  too  easy  to  numerically  integrate  the 
moon’s  heliocentric  coordinates  as  if  it  were  a  planet.  Also,  having  accumulated 
optical  observations  of  the  moon  along  with  those  for  the  planets,  it  seemed  like 
another  small  step  to  utilize  all  these  data  and  differentially  correct  the  lunar  orbit, 
too. 


But  this  move  turned  out  to  be  more  of  a  giant  leap.  The  lunar  orbit 
proved  so  much  more  difficult  to  treat  than  the  planetary  ones  that  we  spent 
almost  all  of  our  time  on  the  moon.  However,  we  will  not  bore  the  reader  with  a 
list  of  difficulties  that  were  encountered. 

The  moon  cannot  be  held  responsible  for  all  the  trouble.  While  there  were  a 
number  of  new  features  to  be  learned,  we  also  made  some  mistakes.  But  none  of 
these  errors  showed  up  in  the  planetary  residuals  and,  hence,  they  were  mostly 
inconsequential  for  that  part  of  the  study.  The  moon  is  much  more  demanding.  Its 
orbit  is  exceedingly  sensitive  to  even  minute  changes  in  the  algorithm  as  well  as  to 
slight  deviations  from  the  proper  procedure  of  running  the  program.  The  latter 
statement  will  be  clarified  later. 

So  all  our  attention  was  focused  on  the  moon  for  a  while.  Fortunately, 
there  were  some  dividends,  predictable  and  unexpected.  First  of  all,  we  think  that 
our  lunar  orbit  is  probably  as  good  as  can  be  obtained  today.  We  feel  strongly 
about  our  approach  to  the  problem.  At  the  present  time,  the  best  and  most 
desirable  procedure  is  a  fit  of  a  precise  numerical  integration  orbit  directly  to 
observations. 
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Secondly,  some  of  the  problems  encountered  prompted  us  to  attempt  an 
improvement  of  the  ephemeris  time  clock,  used  in  the  early  stages  of  this  study. 
Upon  switching  to  atomic  time,  this  program  feature  allowed  us  to  extrapolate  the 
atomic  time  scale  back  to  1912.  Relations  are  given  between  atomic  and  universal  as 
well  as  atomic  and  ephemeris  time.  This  result  alone  should  prove  of  considerable 
interest.  To  the  best  of  our  knowledge,  it  represents  the  first  determination  of  a 
uniform  time  directly  from  observations. 

Finally,  we  obtained  the  coefficient  for  the  tidal  coupling  between  the  orbital 
motion  of  the  moon  and  the  rotation  of  the  earth.  We  find  that  this  phenomenon 
produces  -19"  T2  in  the  moon’s  mean  longitude  with  T  measured  in  centuries. 

v 

A  casual  examination  of  the  problem  may  lead  to  the  conclusion  that  the  six 
lunar  elements,  tidal  coefficient,  and  the  clock  corrections  cannot  all  be  solved  for 
from  a  given  set  of  normal  equations  because  any  signal  in  the  longitude  can  always 
be  removed  by  a  suitable  adjustment  of  the  clock.  This  would  undoubtedly  be  true 
if  we  used  lunar  observations  only  and  imposed  no  other  restrictions.  But  we  also 
have  the  observations  of  sun  and  planets  and  we  have  atomic  time  from  1955  on. 
Since  the  accuracy  'of  atomic  time  far  exceeds  anything  we  could  do  with 

observations,  we  do  not  attempt  clock  corrections  past  1954.  Heuristically,  atomic 

time  from  1955  to  1968  serves  as  the  clock  with  respect  ic  which  the  planetary 

ephemerides  are  constructed.  Prior  to  1955,  the  planetary  observations  strengthen  the 
consistency  of  the  planetary  motions  with  respect  to  one  another  and  with  respect 
to  the  atomic  clock  data  of  1955  to  1968.  The  ephemerides  of  the  inner  planets, 
of  course,  furnish  the  more  sensitive  backward  extensions  of  the  atomic  clock.  The 
lunar  mean  motion  and  change  in  mean  motion  are  fitted  to  the  atomic  clock  and 
to  its  backward  extension  through  the  planetary  ephemerides.  Any  discrepancies  with 
UT,  which  would  be  largest  for  the  lunar  data,  are  then  identified  as  earth  rotation 
inequalities  relative  to  atomic  time. 

We  think  that  the  solution  for  the  planetary  orbits  is  also  a  considerable 

improvement  over  the  classical  ephemerides,  certainly  for  some  planets.  After  all,  the 
residuals  have  been  reduced  to  an  absolute  minimum  in  the  least-squares  sense.  No 
biases  removable  within  the  present  model  remain. 

Among  other  non-orbital  parameters,  we  tried  to  solve  for  the  leading  zonal 
harmonic  of  the  sun.  Results  were  as  follows.  The  formal  standard  deviation  for  J, 
was  1.4  X  10-5,  and  J2  itself  about  the  same  and  with  the  wrong  sign.  Just  about 
the  only  information  which  can  be  extracted  from  that  kind  of  statistics  is  that  J2 
is  not  likely  to  be  much  larger  than  1.4 X  10"5.  A  fair  amount  of  effort  went  into 
this  result.  Since  we  did  not  need  the  relatively  expensive  numerically  integrated 
partial  derivatives,  we  derived  analytical  partials  fiom  Brouwer’s  (1959)  satellite 
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theory.  This  proved  not  as  simple  as  one  might  expect.  We  left  the  fairly  lengthy 
derivation  in  this  paper  for  the  sake  of  documentation.  The  solar  J2  was  a 
suppressed  parameter  in  our  final  iterations. 

As  stated  earlier,  original  plans  called  for  the  addition  of  electronic  data.  They 
were  made  available  to  us  through  the  courtesy  of  Shapiro  and  his  colleagues  at  the 
Lincoln  Laboratory  and  O’Handley  at  the  Jet  Propulsion  Laboratory.  However,  since 
we  had  devoted  so  much  time  to  the  moon,  we  did  not  want  to  delay  the  study 
any  further  and  decided  to  go  for  a  purely  optical  solution;  at  least  for  the  time 
being. 

We  did  not  attempt  to  solve  for  mass  corrections  in  the  present  solution.  There 
were  reasons  for  this  postponement,  over  and  above  the  lack  of  time.  We  had 
already  decided  from  the  beginning  that  numerically  integrated  partial  derivatives 
would  be  used  for  this  purpose.  As  will  be  discussed  in  a  later  section,  inadequate 
partials  in  an  ill-conditioned  normal  matrix  may  lead  to  disaster.  Since  it  is  clear 
that  some  masses  would  be  difficult  to  determine,  the  need  for  precise  differential 
coefficients  was  obvious  and  we  knew  they  would  be  expensive.  However,  as  the 
computer  program  grew,  it  began  taxing  the  storage  capabilities  of  our  computer.  It 
became  inefficient  to  operate,  and  we  priced  our  mass  partials  out  of  reach.  A  new 
program  needs  to  be  written  for  this  purpose. 

Recent  experiments  by  other  researchers  using  various  types  of  partials  tend  to 
support  our  conjecture.  While  there  is  fair  agreement  for  some  masses,  there  is  none 
at  all  for  others.  As  one  would  expect,  the  marginal  masses  of  Pluto  and  Mercury 
are  particularly  troublesome.  In  one  such  set  of  experiments  Pluto’s  mass  fluctuates 
wildly  due  to  only  minor  program  changes.  These  results  strengthen  our  desire  to 
employ  very  accurate  partials.  It  is  clear  that  the  partial  for  the  orbital  elements,  if 
adjured  simultaneously,  must  be  obtained  to  the  same  degree  of  accuracy. 


III.  SPECIAL  FEATURES  OF  DESIGN  AND  ANALYSIS 


The  entire,  program  was  designed  for  and  run  on  our  STRETCH  computer,, 
the  IBM  7030.  The  48  bit  binary  mantissa  yields  14  decimal  places  precision. 
Although  such  a  word  length  is  more  than  adequate  for  most  calculations,  certain 
operations  in  the  integration  routine  are  done  in  double  precision.  They  will  be 
pointed  out  later.  Overall  design  precision  is  0"01.  Practically  all  calculations  meet 
this  criterion.  However,  the  truncation  error  in  the  lunar  orbit  slightly  exceeds  this 
threshold.  As  will  be  explained  later,  the  maximum  effective  truncation  error  here  is 
0"013. 

The  starting  time  of  the  integration  is  JD  242  0000.5  which  corresponds  to 
1913  August  21.0.  The  reasons  for  this  choice  are  given  in  the  section  on 
observational  material.  Observations  span  the  20,000  days  to  JD  244  0000.5.  In  case 
of  the  moon,  the  data  were  separated  into  two  groups,  those  before  and  after  1925. 
As  it  were,  all  lunar  Six-Inch  observations  have  recently  been  uniformly  reduced  by 
Adams,  Klock,  and  Scott  (1969).  They  were  put  on  the  FK4  system  and  have  all 

received  limb  corrections.  Since  no  such  uniform  treatment  exists  for  the  data 

before  1925,  the  two  sets  should  receive  different  weights.  The  weights  calculated 
and  employed  are  found  in  Table  IX. 

Among  the  various  perturbing  forces  studied,  but  not  discussed  elsewhere  in 
this  paper,  is  the  solar  radiation  pressure.  Although  small,  this  acceleration  is  known 
to  cause  problems  in  the  case  of  some  artificial  satellites.  In  our  problem,  the 

effects  are  entirely  negligible.  Peters  (1964)  has  pointed  out  that  radiation  pressure 
on  any  given  planet  acts  so  as  to  decrease  the  apparent  mass.  Even  for  Mercury,  for 
which  the  effect  would  be  largest,  solar  radiation  pressure  corresponds  to  a  change 
of  its  mass  by  only  2X  10'7. 

The  partial  derivatives  used  in  our  program  were  obtained  partly  by 

numerical  integration  but  mostly  analytically.  They  are  discussed  in  later  sections. 
For  some  experiments  we  also  approximated  partials  of  the  form  3D/9P  by  finite 
difference  ratios  AD/AP. 

We  made  some  effort  to  include  the  lunar  tesseral  harmonics  of  degree  and 
order  two  in  our  model.  All  runs  with  such  terms  resulted  in  some  peculiar  signals. 
Since  we  were  unable  to  find  an  error  in  either  formulation  or  computer  program, 
these  terms  were  eventually  suppressed  by  putting  the  coefficients  c22  and  s22  to 
zero.  The  derivation  of  the  corresponding  accelerations  is  retained  in  this  report. 
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It  may  seem  that  the  documentation  of  derivations  is  not  at  all  uniformly 
complete.  There  is  a  reason  for  this.  We  saw  no  need  to  record  anything  that  is 
well  known  and  easily  found.  We  did,  however,  try  to  retain  the  essential  steps  of 
new  developments  or  those  derivations  which  we  always  have  difficulties  finding. 

The  reader  will  see  in  a  later  part  of  this  paper  that  not  all  73  parameters 
can  be  solved  for  simultaneously.  They  are  divided  into  two  subsets,  one  containing 
the  planetary  elements  and  the  other  the  lunar  elements  and  associated  parameters. 
These  sets  are  alternately  solved  for  in  consecutive  least-squares  solutions.  The  formal 
standard  deviations  produced  with  each  such  partial  solution  are  ignored.  After  final 
convergence  we  made  a  global  solution  only  in  order  to  obtain  a  more  realistic 
variance-covariance  matrix.  All  standard  deviations  are  taken  from  this  matrix. 

The  same  matrix  was  used  to  get  the  unitized  correlation  coefficients.  Some 
of  these  coefficients  are  tabulated  near  the  end  of  the  paper.  It  is  clear  that  only 
limited  information  can  be  obtained  from  such  figures. 


IV.  OBSERVATIONS 


As  remarked  earlier,  we  had  plans  of  collecting  observations  of  sun,  moon, 
and  planets  back  to  1750  or  so.  A  preliminary  library  search  revealed  that  this  task 
was  probably  even  more  formidable  than  expected.  The  major  difficulty  was  due  to 
the  fact  that  in  older  publications  the  observations  of  the  bodies  of  interest  to  us 
were  listed  among  those  taken  on  stars.  Starting  about  1911,  though,  results  of  the 
U.  S.  Naval  Observatory  transit  circle  programs  were  published  separately  for  each 
planet.  It  seemed  to  us  that  the  data  from  1911  on  were  a  good  starting  set  for 
our  work.  Shapiro  (1970)  recently  pointed  out  that  the  USNO  transit  circle  results 
are  of  higher  quality  than  any  other  series. 

In  the  meantime,  Shapiro  has  actually  gone  back  about  200  years  and 
collected  meridan  circle  observatories  of  the  solar  system  bodies  taken  at  several 
major  observatories.  The  product  of  this  tremendous  and  commendable  effort  are 
something  like  350,000  or  400,000  observations  in  right  ascension  and  declination. 

Parts  of  the  initial  data  set  punched  by  us  were  eventually  replaced  with 
updated  reductions  done  at  the  USNO.  However,  we  presently  plan  to  retain  all  of 
the  observational  material  in  machine  readable  form  for  a  few  years.  Hence,  the 
data  as  originally  published  as  well  as  those  eventually  used  by  us  in  our  orbit 
improvements  will  be  available  to  other  interested  researchers. 

The  original  observations  were  found  in  five  volumes  of  the  Publications  of 
the  United  States  Naval  Observatory,  Second  Series  (1927,  1933,  1948,  1949,  1952). 
Table  VII  lists  the  volumes  we  punched.  All  these  data  were  taken  with  the 
Six-Inch  or  Nine-Inch  transit  circle.  The  “source  code”,  given  in  the  last  column  of 
Table  VII,  is  a  ready  identification  of  the  particular  observation  series,  and  it  is 
included  on  every  observation  card  punched. 

Table  VIII  represents  an  attempt  of  listing  the  sources  of  all  the  observations 
used  in  our  final  solution.  As  may  be  seen,  three  agencies  cooperated  in  punching 
the  data  onto  cards.  Sections  of  our  master  tape  were  updated  a  number  of  times 
when  “better”  observation,  usually  the  results  of  a  more  definitive  reduction,  became 
available.  Such  changes  were  made  for  all  planets  in  a  given  time  span  or,  as  in  the 
case  of  V  XIX  Pt  III,  for  one  body  at  a  time.  After  so  many  revisions,  minor 
discrepancies  between  the  publications  and  our  master  tape  are  likely.  The  only 
major  discrepancy  we  are  aware  of  is  an  apparent  loss  of  340  days  of  data  from 
Circular  No.  1 1 5(C1 1 5).  It  is  not  known  how  this  occured  and  no  attempt  will  be 
made  at  this  time  to  recover  the  missing  points.  No  data  was  lost  from  Cl 27  as 
the  cutoff  was  planned. 
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All  observation  cards  are  now  available  in  one  unique  format.  This  was  not 
so  at  first.  In  punching  from  the  USNO  publications,  some  15  to  20  differen 
formats  were  encountered.  It  was  clear  that  all  cards  had  to  be  converted  to  a 
uniform  format.  In  checking  with  JPL  in  the  fall  of  1967,  we  found  that  they  had 
just  devised  a  card  format  for  in-house  usage.  With  a  few  minor  modifications,  the 
format  was  made  to  fit  the  needs  of  JPL,  NWL,  and  USNO.  Other  interested  partie. 
were  asked  for  inputs,  and  a  standard  format  emerged.  The  latter  is  explained  in 
detail  by  O’Handley  (1968). 

A  few  additional  comments  have  to  be  made  about  Pluto.  Since  this  object 
is  too  faint  to  be  observed  with  a  transit  instrument,  all  observations  are 
photographic.  The  entire  material  used  in  our  Pluto  orbit  correction  (Cohen, 
Hubbard,  Oesterwinter,  1967)  was  employed  here  again.  However,  all 
Sharaf-Budnikova  normal  places  were  treated  as  individual  observations.  We  believe 
that  this  simplification  is  of  little  consequence  in  the  present  investigation.  Added  to 
this  data  set  were  two  Russian  observations  taken  in  1965  (Chernykh  and  Chernykh, 
1967).  The  source  code  D  XX  for  Pluto  in  Table  VIII  is  only  the  general  form  of 
seven  different  codes  used.  They  are  explained  in  Table  IX. 

Our  master  tape  also  contains  all  the  observations  made  on  Ceres,  Pallas, 
Juno,  and  Vesta  as  listed  in  D7,  D9,  Dl,  D2,  D3  and  D4.  These  data  were  not 
used. 


Observation  times  given  in  any  of  the  publications  were  ignored.  We 
calculated  our  own  observation  times  from  the  obseived  right  ascensions  using  a 
procedure  described  below.  Hence,  any  observation  of  declination  alone  or 
incomplete  in  right  ascension  is  useless  and  was  disregarded.  On  the  other  hand,  the 
master  tape  contains  observations  with  missing  or  incomplete,  information  in 
declination.  Although  still  useful  in  right  ascension,  such  data  points  were  eliminated 
at  a  later  stage  in  an  effort  to  simplify  the  program. 

Most  observations  showing  large  residuals  on  preliminary  nins  were  also 
removed  from  the  master  mpe.  In  view  of  the  large  amount  of  data  available  no 
effort  was  made  to  trace  the  discrepancies.  Simple  keypunch  errors  are  the  likely 
cause  in  many  cases. 

First  runs  showed  a  distinct  bias  in  the  right  ascension  residuals  for  early 
observations  of  sun,  Mercury,  and  Venus.  The  problem  was  quickly  traced  to  a 
difference  in  procedure  employed  in  -the  reduction  of  source  code  6  data  in  contrast 
to  later  volumes.  The  so-called  equinox  correction  of  - 1  "21 8  was  not  applied  to  the 
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observed  right  ascension  in  those  days.  We  simply  subtracted  0?  08 1  from  all  right 
ascensions  in  source  code  6.  The  data  on  our  present  master  tape  are  thus 
corrected. 

This  modification  was  the  only  change  made  in  any  of  the  original 
observations.  Since  the  master  tape  is  saved,  the  data  are  preserved  in  this  form. 

There  are  a  number  of  other  corrections  that  have  to  be  applied  but  none 

of  these  have  been  added  to  the  original  master  tape.  The  old  astronomical  tradition 

of  leaving  observed  places  alone  is  a  very  sensible  one.  The  exception  discussed 

above,  however,  is  where  it  belongs,  we  think,  since  it  removes  an  inconsistency  in 

the  original  data  reduction. 

The  reader  may  wish  to  check  the  section  on  residuals  for  additional 
corrections.  Among  others,  we  made  a  bias  correction  just  like  the  one  above.  The 
important  difference  is  that  we  were  unable  to  find  a  cause  for  the  discrepancies. 
Hence,  we  did  not  wish  to  alter  the  original  data. 


V. 


DATA  PREPARATION 


The  observations  on  the  master  tape  cannot  yet  be  compared  with  computed 
places.  A  number  of  modifications  are  still  to  be  made.  Unless  otherwise  stated,  the 
following  discussions  do  not  pertain  to  Pluto  observations. 

Older  catalogs  did  not  use  universal  time.  Since  approximate  times  are 
needed,  as  will  be  seen  in  a  moment,  we  first  converted  all  published  observations 
times  to  UT.  Table  I  shows  the  scheme.  The  reader  will  recognize  immediately  that 
0?5  of  these  corrections  is  due  to  the  fact  that  the  astronomical  day  before  1925 
began  at  noon.  The  value  0?2  is  simply  the  reduction  from  EST  to  UT.  These 
numbers  are  required  only  to  the  nearest  0?1. 


TABLE  I 

Correcting  Catalog  Times  to  UT 


Source  Code 

Add  to  t 

6 

0?7 

7(t  <  1924) 

0.7 

7 (t  >  1925) 

0.2 

8,  9,  0 

0.2 

The  derivation  of  observation  time  from  the  recorded  right  ascension, 
subsequently  labelled  a,  is  of  utmost  importance  to  the  success  of  the  entire 
operation.  Hence,  we  will  record  here  the  essential  steps. 

Consider  any  calendar  date.  Let  I  represent  the  year,  J  the  month,  K  the 
integral  part  of  the  day,  and  d  its  fractional  part.  Then  the  Julian  Date  at  OnUT 
can  be  obtained  from: 


JD0  =  1721074 +  K  + 1461  *(1  +  (J-  14)  /  12)  /  4  +  367  *  (J  -  2  -  12* 
((J  -  14)/  12))/  12 +  (24002-  12  *  I  -  J)  /  1200  -  0.5 . 


(1) 
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Except  for  the  last  term,  0.5,  this  is  obviously  a  Fortran  statement  employing 
integer  arithmetic.  It  can  be  used  at  once  on  almost  any  computer.  This  formula  is 
based  on  one  originally  supplied  by  Seidelmann  (1967)  at  the  USNO.  We  modified 
the  relation  so  that  it  is  now  valid  in  the  open  interval  1600  to  2100. 

Next  we  need  Greenwich  apparent  sidereal  time  at  0hUT: 


GAST0  =  05*2764  9045  +0.002  737  909  298  (JDQ  -243  1090.5) 

+  0?000  000  7716  A*  cose. 


(2) 


This  number  must  be  between  0  and  1  which  can  be  done  by  adding  and 
subtracting  multiples  of  1.  The  last  term  is  usually  called  equation  of  the  equinoxes. 
Here  ASk  is  the  nutation  in  longitude,  expressed  in  seconds  of  arc.  As  in  other 
places,  e  is  the  obliquity  of  the  ecliptic. 

With  this  precomputation,  the  fractional  part  of  the  day  follows  from 


fd  =  0.997  269  566  (ad  +0d2140  7233  -  GAST0) .  (3) 


If  fd  does  not  fall  into  the  interval  from  0  to  1,  the  expression  in  parentheses  is 
increased  or  decreased  by  one.  In  this  formula,  a  is  the  observed  right  ascension, 
and  the  next  term  represents  the  longitude  of  the  USNO.  Finally,  the  observation 
time  in  UT  is 


tUT  “  JDQ  +  fd 


(4) 


The  procedure  does  not  always  yield  a  unique  time  inasmuch  as  two 
observation  times,  about  24h  apart,  are  possible.  The  ambiguity  is  easily  resolved  as 
long  as  the  original  observation  time  is  recorded  to  the  nearest  0?  1  or  so. 
Comparing  tUT  wiih  the  approximate  JD0+d  eliminates  the  spurious  solution. 
However,  in  source  code  6  times  are  given  to  the  nearest  day  only.  In  this  case  we 
wrote  both  possible  tUT  on  the  master  tape  and  rejected  the  erroneous  one  only 
after  inspection  of  the  residuals. 
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It  should  be  emphasized  that  we  obtained  observations  times  for  all  data, 
except  Pluto,  in  the  manner  described  above.  In  other  words,  this  applies  not  only 
to  data  punched  here,  but  also  to  the  more  recent  data  sets  provided  by  JPL  and 
USNO. 


At  this  stage,  only  the  Pluto  observations  are  still  in  calendar  time. 
Conversion  to  the  Julian  data  is  accomplished  using  (1)  and,  since  the  exact 
fractional  day  is  available, 


tUT  =  JD0+d.  (5) 

The  discussion  up  to  this  point  explains  the  data  written  on  the  “master 

tape”.  The  latter  is  now  fed  into  our  data  conversion  program  (DCP)  which 

produces  the  “observation  tape”  actually  used  by  the  main  program. 

All  transit  circle  observations  are  referred  to  the  true  equator  and  equinox  of 
date.  Since  the  numerical  integration  yields  results  in  a  fixed  frame,  the  observations 
must  be  corrected  for  nutation  and  precession.  We  selected  the  mean  frame  of 
1950.0  to  which  the  Pluto  data  are  already  referred. 

Most  of  the  algorithm  for  calculating  nutation  and  precession  was  lifted 
directly  from  the  Explanatory  Supplement  (1961).  The  required  elements  of  sun  and 
moon  are  given  there  on  p.  44.  The  trigonometric  series  for  A'P  and  Ae,  nutation 
in  longitude  and  obliquity,  pp.  44  and  45,  were  truncated  after  terms  with 

coefficient  0"0050.  We  also  copied  the  equation  for  e,  the  mean  obliquity  of  the 
eccliptic  of  date,  written  as  a  polynomical  in  time,  from  p.  98.  With  these 
preparations,  we  calculate  the  effects  of  nutation  as  follows.  Adding  Ae  to  e  gives 
the  true  obliquity  of  date.  Together  with  the  available  a  and  o,  we  now  obtain  the 
true  longitude  and  latitude  from  the  spherical  trigonometry  relating  equator  and 

ecliptic.  It  is  only  necessary  to  subtract  A'P  from  the  true  longitude,  and  one  has 
mean  longitude  and  latitude  of  date.  Together  with  e,  this  time,  one  solves  for  aM 
and  SM,  where  the  subscript  M  indicates  values  referred  to  the  mean  frame  of  date. 
In  order  to  be  generally  applicable,  such  angles  arc  always  obtained  from  inverse 
tangent  relations. 

Also  copied  from  the  Explanatory  Supplement  are  the  expressions  for 
precession.  The  quantities  f0,  z,  0  are  given  on  p.  30.  We  rewrote  the  spherical 
relations  on  p.  31  a  bit.  Evidently  one  can  combine  two  equations  to  obtain  tan 
(a-z),  solve  for  (a-z),  and  then  for  a.  Subsequently  one  can  also  get  6  from  an 
arc  tan. 
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All  times  are  still  referred  to  the  UT  clock.  The  last  function  of  our  DCP  is 
the  conversion  of  the  observation  times  to  the  adopted  clock.  For  this,  ephemeris 
time  was  chosen  in  our  earlier  solutioi  .  The  AT’s  to  be  added  to  tUT  were  then 
obtained  from  the  American  Ephemeris  and  from  Brouwer’s  (1952)  original  paper  on 
the  subject.  Later  we  switched  to  atomic  time.  Here  the  quantity  added  from 
1955.5  on  was  32?15  +  A.l  -  U.T.2,  given  in  the  American  Ephemeris.  All  values 
before  1955.5  were  obtained  as  part  of  our  solution.  More  about  this  will  be  said 
in  a  later  section. 
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VI.  ALGORITHM  AND  DERIVATIONS 


1.  PRECOMPUTATIONS 

One  feature  of  the  program,  exercised  only  in  special  tests,  permits  to 
designate  which  bodies  gravitationally  affect  any  given  planet.  This  option  could  be 
used  not  only  for  very  small  asteriods  or  artificial  bodies,  but  also  in  cases  where, 
for  example,  Pluto’s  effect  on  Mercury  is  negligible.  Using  the  list  of  perturbing 
planets,  separately  for  each  body  being  integrated,  the  program  selects  the 
appropriate  equation  of  motion,  with  the  minimum  number  of  terms  required,  from 
available  series  of  such  equations. 

In  case  the  initial  conditions  of  any  planet  are  referred  to  the  mean 
equator  and  equinox  of  T  different  from  the  selected  T0 ,  the  updating  is  quickly 
made.  The  angles  £0  ,  z,  and  0  are  obtained  as  discussed  above.  Then  the  precession 
matrix 


(cos  f0  cos  6  cos  z  -  sin  f0  sin  z 
cos  f0  cos  6  sin  z  +  sin  cos  z 
cos  sin  0 


-  sin  f  0  cos  0  cos  z  -  cos  f0  sin  z 

-  sin  f0  cos  0  sin  z  +  cos  f0  cos  z 

-  sin  f0  sin  0 


sin  0  cos  z 
sin  0  sin  z 
cos  0 


Finally 


?(r0)  =pf(x) 
and 

?(T0)  =PF(t) 

Note  that  T0  =  JD  2433282.423  =  B.Y.  19500  in  all  our  production  work. 

The  obliquity  of  the  ecliptic  is  required  in  various  places.  The  mean 
value  if  given  by 

e  =  23?452294-  0?0130I25AT- 0?00000164AT2 +0?0000005030  AT3 
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where 


AT  = 


T0  -  2415020.0 
36525 


2.  NUMERICAL  INTEGRATION  ROUTINE 

The  differential  equations  described  in  the  next  section  were 
numerically  integrated  using  our  present  “Cowell”  routine.  This  program  is  the  latest 
modification  of  a  series  of  such  integrators  developed  and  improved  over  the  years 
by  Cohen  and  Hubbard.  Some  versions  are  described  by  Hubbard  and  Broadwater  in 
the  -  open  literature  and  in  various  local  reports.  We  used  them  in  several  studies  of 
planetary  motion,  and  they  are  being  employed  extensively  in  satellite  geodesy. 

The  program  begins  by  deriving  a  number  of  auxiliary  coefficients 
using  the  Cowell  and  Adams  coefficients  Cj  and  a;  given  in  Table  X. 


TH)i( J) 


j  =  0,l,2,“*,a  +  b 


where  [j  +  l,b]  is  the  smaller  of  j  + 1  and  b, 


■  I  a  +  b  +  1 
b:  =  (-DM 
J  \  j+i 


t  ■  H)J40a' 


i  =  0,l,2,**',a  +  b 


j  -  0,l,2,,,*,a  +  b 


j  =  0,l,2,***,c 


j  =  0,1, 2, -".c 


cf  =  '  2  c. 
1  j-o  1 


i  =  0,l,2***,m 


j  =  0,l,2,***,m 


f 


i  =  0,1 ,2,*  •  *,m 


i 


a!  =  2  a. 

1  j=o  J 

*  ■  (-'MC> 


j  =  0,l,2,***,m 

j  =  0,1 ,2,*  •  *  ,m 

j  =  0,1,2,*  **,m 


Here  c  and  m  are  the  order  of  starting  and  running  routine,  respectively.  The 
controls  a  and  b  are  better  explained  farther  below. 

The  starting  sequence  begins  by  evaluating  the  differential  equations  at 
epoch,  that  is,  r Q.  The  starting  table  is  initialized  by  putting 


rn  =  rQ  for  -a  <  n  <  b 


Figure  1  depicts  the  scheme.  In  words,  the  unknown  accelerations  are  set  equal  to 
the  initial  value  in  the  range  from  -a  to  b.  The  numbers  a  and  b  cannot  exceed  c. 
Next  the  position  on  the  time  line  before  epoch  is  given  by 

.  a+b 

r-l  =  VhF0+h  i=£0Ti  fb-i  (6) 


where  h  is  the  integration  step  size. 


The  accelerations  r  are  extrapolated  beyond  -a  and  b  to  c  lines  before  and  after 
epoch  by 


a+b 

i?0  bi  *-*'*J 


r_  = 


a+b 

2  r„  .  j 
j=  0  J  n“1_J 


n  =  a  +  l,a  +  2,***,c 
n  =  b  +  l,b  +  2,***,c 
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FIGURE  1 

Scheme  for  Generating  Starting  Table 
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Position  and  velocity,  the  latter  always  required  in  our  model,  are  made  in  the 
range  a  to  b  only: 


=  2r„  ,  -  rn  ,  +  h2  Z1  a,*  rn  • 
n-l  n-2  j_0  j  n-j 


n  =  1,2 


>b 


F_n  =  2r.n+1  -  r.n+2+h2.Sa*f_n+j 


r  - 


"A1*)  f-"*‘ 


n  =  2,3,**,,a 
n  =  1 ,2, *  *  *,b 

n  =  l,2,*'#,a 


Finally,  the  accelerations  in  the  same  interval  are  calculated  by  evaluating  the 
differential  equations  using  the  coordinates  just  generated. 

At  this  point  the  program  tests  whether  K  cycles  through  the  starting 
routine  have  been  completed.  We  used  K  =  4  for  almost  all  of  our  runs  based  on 
previous  experience.  If  K  is  not  yet  reached,  the  computer  returns  to  (6)  for  a  new 
cycle.  When  K  cycles  are  finished,  a  convergence  test  is  made.  For  this  purpose,  the 
absolute  difference  in  the  accelerations  on  two  successive  sweeps  is  monitored  for 
each  differential  equation.  When  the  maximum  of  these  numbers  reaches  a  minimum, 
the  starting  tables  are  considered  converged. 

The  running  procedure,  of  course,  advances  one  step  at  a  time. 
Coordinates  on  the  next  time  line  are  made  by  the  predictor  formulae 


=  2  rn  ,  -  fn  ,  +  h2  2  <*:?.; 
n-l  n-2  j  n-i-j 


m 


r „  =  r„  .  +  h2  ft  r„  ,  ■ 

n  n-l  j=o  J  n-1_j 


where  m  is  the  order,  in  general  different  from  c.  The  corresponding  accelerations 
arc  calculated  from  the  differential  equations.  The  routine  may  be  exited  at  this 
point  or  coordinates  can  be  refined  with  the  corrector  formulae 
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Again,  new  accelerations  are  computed.  In  order  to  enhance  stability,  the  last  step  is 
always  the  evaluation  of  the  accelerations. 

After  extensive  testing,  we  used  the  following  constants  for  all  our 
production  work: 


h  =  0?4 
c  =  14 
a  =  b  =  7 
m  =  12 


For  the  planets  Mars  through  Pluto,  it  was  found  that  above  predictor  formulae 
sufficed.  For  all  other  differential  equations,  including  the  various  variational 
equations,  the  predictor-corrector  loop  was  required. 

In  order  to  minimize  roundoff,  coordinates  were  always  computed 
and  stored  as  double  precision  numbers.  We  examined  this  device  quite  recently  and 
found  that  in  many  applications  single  precision  calculations  will  give  identical 
results. 

3.  EQUATIONS  OF  MOTION 

There  is  no  need  to  reproduce  the  particle  terms  of  the  accelerations. 
As  indicated  before,  the  moon  is  treated  like  any  other  planet,  that  is,  it  has  its 
own  heliocentric  orbit.  This  choice  makes  for  a  simpler  program,  and  no  numerical 
disadvantages  result. 

Evaluating  the  differential  equations  is  usually  the  most 
time-consuming  operation  in  a  numerical  integration  routine  such  as  this.  Hence,  we 
took  pains  to  rewrite  the  principal  equations  of  motion  in  a  form  advantageous  to 
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the  computer.  Also,  if  a  body  of  negligible  mass  was  integrated,  its  effects  on  the 
other  planets  was  not  simply  suppressed  by  putting  the  mass  to  zero,  but  the  terms 
containing  its  mass  were  altogether  absent. 

Relativistic  terms  could  be  computed  or  by-passed  separately  for  each 
planet  through  a  simple  input  control.  Except  for  experiments,  they  were  always 
calculated  for  every  planet.  The  equations,  as  quoted  by  Brouwer  and  Clemence 
(1961),  are  based  on  the  standard  Schwarzschild  metric.  Experts  assuie  us  that  these 
approximations  are  more  than  adequate  to  process  optical  observations,  as  Brouwear 
and  Clemence  stated.  Rewritten  a  bit,  the  additional  terms  due  to  general  relatively 
are 


2p  r  t  ^  •  2  (  c2r  \ 

Ar  = - !  j/r  -  r2v2  +  (r  •  P)2  (1  + - \  r 

c2!-5  |  \  2(c2r-2ju)/ 


where 


2 fi  ?  *  r  i 
+  —  - r 

r2  c2r-2M 


H  =  k2(l  +m) 


and 


v  =  IF  1 


Spherical  harmonics  of  earth  and  moon 

It  is  clear  that  some  of  the  deviations  of  earth  and  moon  from 
spherical  symmetry  play  a  role  in  our  work.  The  effects  of  zonal  harmonics  can 
easily  be  assessed  from  existing  artificial  satellite  theory,  such  as  Brouwer’s  (1959). 
For  some  of  the  other  terms  we  made  a  series  of  computer  runs,  with  and  without 
the  pertinent  accelerations,  and  compared  results. 

The  forces  due  to  the  various  gravitational  terms  were  derived  as 
follows.  In  inertia]  space,  the  force  potential  F  for  a  spherical  sun  and  nonspherical 
earth  and  moon  can  be  written  as 
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(8) 


rrij  m2 


L»  "  ^12(^1^12)“  ^21  (^2^1 2^) 


momi 


“  [l-U10(Ppf0i)]+T^  [l-U20(i'2,?02)]+0(UX 

*01  r02 


U) 


The  subscripts  0,  1,  and  2  refer  to  sun,  earth,  and  moon.  The  is  the  perturbing 
potential  of  body  i  affecting  body  j,  and  Pj  is  a  vector  along  the  spin  axis  of  body 
i.  The  meaning  of  the  other  symbols  seems  clear. 

There  is  no  need  to  consider  the  oblate  figure  of  the  sun  in  this 
development.  It  will  be  treated  when  required  in  a  later  section.  In  fact,  if  the  sun 
were  left  out  altogether,  the  resulting  equations  would  still  be  found  to  provide  very 
good  approximations. 

We  will  employ  the  usual  development  of  (8)  using  essentially  the 
notation  adopted  by  the  IAU  (Hagihara,  1962).  However,  we  shall  dispense  with  the 
commas  in  subscripts.  Moreover,  gravitational  coefficients  belonging  to  the  earth  will 
be  denoted  by  capital  letters,  those  of  the  moon  by  lower  case  letters.  Equation  (8) 
then  becomes 


1  m.m,  f  «°  n  /RA"  _ 

=  ~  LI+  n=.  4^)Pn(Sin',2)X(CnmC0SmXlJ+SnmSinmXl2) 

+  nJ,  m?o(r^j P"  (Si"  ^l)X(Cnm  cos  m\  l  +  Snmsin  mX2 ]  >] 


(9) 


mnm.  r  »  n  /R,\n  „ 

—  1+n?!m?0rp>n^o)x(Cn 

roi  L  1  m_0  \r0 1/ 


„  cosmXin  +  S _ sin  mX 

m  10  nm 


10^ 


■I  in0m2  ['  ~  n  /RA" 

'o,  hn^?o(3C(Sin^)X(C" 


nicosmX20  +  snm  sinmX20) 


) 
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In  this  formula,  R  stands  for  the  equatorial  radius  and  P  for  the  associated 
Lengendre  polynomials.  The  symbols  0-  and  XSj  are  the  latitude  and  longitude  of 
body  j  in  an  equatorial  frame  of  body  i.  Since  the  origin  of  the  customary 
rectangular  frame  is,  hopefully,  at  the  center  of  mass,  terms  with  n  =  1  vanish. 
Also,  C21  =  S21  =0  and  c21  =  s21  =0  since  the  z-axis  is  taken  along  the 
principal  axis  of  inertia.  From  various  considerations  we  know  that  we  need  the 
terms  with  C20and  C30,  but  those  factors  by  C22  and  S22  are  negligible.  In  case 
of  the  moon,  we  need  to  go  to  degree  two  only,  but  c22  and  s22  wiil  be  retained. 
Recall  that  the  earth  is  a  synchronous  satellite  over  the.  first  meridian  of  the  moon. 
Finally,  the  central  terms  are  already  in  the  equations  of  motion  so  that  only  the 
distrubing  part  AF  of  (9)  is  required.  There  remains 


1 

rj  AF  =  m!m2 


C 


+  m,m2  ~ 

12 


P?(sin0. 
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)XC20  +  P2^S‘n<32l)X 


c22  cos  2X2  j  +s22  sin  2X2 


.) 


(10) 


+  nlom,  XC20  +  ^-  P°(sin  01O)  X  C  J 

01  L  oi  J 

+  m0,n2  'T'[P?(Si"  X  c20  +  P2(Sin  X  (C22  C0S  2X20  +  S22  Si"  2X2ol 
02  L  J 


The  equations  of  motion  in  inertial  space  would  be  of  the  form 


'Vi 


=  Vi 


i  =  0,1,2 


where  the  tilde  indicates  differentation  with  respect  to  inertia!  coordinates.  Since  we 
want  the  heliocentric  coordinates 


the  equations  of  motion  in  the  latter  frame  are  found  to  be 


j  =  1,2 
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-  Ur  +rIVF+r^ 


m,  ..m0 


m„ 


/  1  l  \  1 

Ar,  =/—  +—  V.AF+ —  V.AF 

\m2  m0/  m0 


(10 


The  gradient  now  indicates  differentation  w.r.t.  heliocentric  coordinates. 


Execution  of  (11)  is  a  bit  lengthy  but  fairly  straight-forward.  In  fact, 
Vj  AF  is  available  from  many  different  sources.  In  doing  V2  AF,  the  rapid 
precessional  rate  of  the  lunar  equator  has  to  be  taken  into  account,  and  the 
transformation  into  the  adopted  frame,  the  earth’s  equator,  has  to  be  made.  There 
are  only  a  few  details  in  the  development  we  wish  to  record  here. 

Now  that  we  are  in  the  heliocentric  frame,  we  will  write 


forr, 


01 


and 


r2  for  F0? 

We  will  also  use  the  geocentric  vector  to  the  moon 


r 


G 


It  is  then  easily  seen  that 


and 


■  „  zg 

sinpt2  =  — 


sin/3 
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Both  these  relations  are  needed  in  (10). 


Let  us  now  turn  our  attention  to  the  task  of  relating  coordinates 
referred  to  the  equator  of  the  moon  to  that  of  the  earth.  For  a  moment,  let  any 
vector  without  subscript  be  presented  in  the  latter  frame,  then  the  familiar  rotation 


rECL 


0  0 
cose  sine  lr 
-sine  cosey 


(12) 


gives  the  same  vector  in  the  frame  of  the  ecliptic.  We  will  use  e  =  23?446.  With 
the  aid  of  Figure  2,  in  the  frame  of  the  lunar  equator  this  vector  becomes 


rM 


'1  0  0  \  /cosH 

=  I  0  cos  I  sin  I  1  /-  sin  H 

,0  -sin  I  cos  I  A  0 


sin  H 
cos  H 
0 


(13) 


Here  I  is  the  inclination  of  the  lunar  equator  on  the  ecliptic  for  which  we  take  the 
mean  value 


I  =  1?535. 


H  stands  for  the  ascending  node  of  the  lunar  equator  on  the  ecliptic.  Again,  for 
our  purpose  we  need  not  be  concerned  with  periodic  perturbations.  Hence,  we  find, 
appealing  to  lunar  theory  and  Cassinis  laws, 


H  =  79?  183275-  0?05295  39222  (t  -  241  5020.0) 


where  t  is  expressed  as  a  Julian  date.  Putting  (12)  into  (13),  wc  have  a  relation  of 
the  form 


rM  =  M  r. 


(14) 
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LUNAR  EQUATOR 


FIGURE  2 


Relation  of  Lunar  Equator  to  Ecliptic 


Upon  substitution  of  the  numerical  values  given,  one  finds 


(cos  H 
-.99964  sin  H 
0.02679  sin  H 


0.91744  sin  H 
0.91711  cos  H -.01066 
-.02458  cos  H -.39774 


0.39788  sin  H 
0.39774  cos  H  +.02458 
-.01066  cos  H +.91711 


(15) 


Now  we  are  ready  to  consider 


and 


sin  02 , 


sin  02O 


(16) 


Because  of  (14),  we  can  write 


?2i  =  M(r,  -  r2)  =  -MrG  (17) 

r20  =  M(0-  r2)  =  -Mp2 

Let  now  the  last  row  of  M  be  designated  M3,  that  is 


(18) 


Then,  because  of  (17), 
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and 


Z20  ~  "M3  *  f2 


Putting  this  into  (16),  we  have  the  desired  result 


sin  P2 1 


sin  p2o  = 


In  order  to  relate  the  two  selenocentric  longitudes  in  (10)  to  available 
coordinates,  we  begin  with  Figure  3.  The  unprimed  coordinates  are  already  familiar. 
The  primed  ones  refer  to  a  frame  rotating  with  the  moon  oriented  such  that 
xM  goes  through  the  first  meridian  in  the  equator  from  the  latter’s  ascending 
node  on  the  ecliptic.  The  angle  X  is  the  selenocentric  longitude  of  an  object  at  S. 
The  relation  between  any  object’s  polar  and  rectangular  coordinates  are  obviously 


rcos/3  cos  (X  +0)' 
r|  cos/3  sin  (X  +  0) 
sin/3 


Hence 


X 
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-  0 


X20  =  tan' 


-l 


y  20 
*20 


-  0 


(19). 


Now  consider  (17).  If  we  designate  the  first  two  rows  of  M  by  M,  and  M2, 
analagous  to  (18),  there  follows 
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FIGURE  3 


Relation  of  Selenocentric  Longitude  to  Ecliptic 


x21  -Mj  •  rG 

*21 

=  -M. 

X2  0  =  ~^1  *  *2 

*  20 

=  -M. 

With  these  relations  (19)  becomes 


X 


21 


tan-1 


-_M,  *  fr. 

-Mj  -rG 


^20 


All  quantities  in  (10)  are  now  expressed  in  terms  of  available 
coordinates  so  that  the  differentiations  indicated  in  (11)  can  be  performed.  The  final 
results  together  with  an  outline  of  the  sequence  of  calculations  follows. 

Notation  is  often  a  problem,  and  this  case  is  no  exception.  While  it  is 
convenient  in  this  section  to  let  the  subscript  1  refer  to  the  earth  and  2  to  the 
moon,  the  subscripts  actually  used  in  the  computer  program  were  3  and  11.  Since 
it  will  be  necessary  to  use  the  latter  elsewhere  in  this  paper,  this  minor 
inconvenience  should  be  noted. 

Calculations  begin  with  (14)  followed  by  the  various  in  (15).  Next 

we  find 


e  =  0o+w(t- 10) 


Here  <o  is  the  draconic  rate  of  the  moon.  The  numerical  values  employed  by  us 
were 


cu  =  13? 2293  5045 /day 


31 


/ 


and 


0O  =  25?531 


The  longitudes  X21  and  X20  follow  from  (20).  Lef 

'O' 

0 
1 


.  p=l 

and  calculate 

Si  ="2k2c20Ri 

ii 

ZG 

r2  “  1 
lG 

)  Fo  ~  3zg 

S2  =  2k2C20R? 

?R 

72 

5  ~2~  ~  1 
r?  , 

)r|-2z,f] 

n  1  7 

s;  - 


C30ri  rs  [5rf  (7  ^  '  fr,  -  3 (5 ri  - 


^=;*2c3.Rf^f 


(c22  COS  2?lj 
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Finally,  the  additions  to  the  accelerations  of  earth  and  moon  become 

Ari  =  m2  (si  +  S3  +  S5  ”  S6  “  S7  +  S8  “  S9  +  S1  o)  +  0  +  ml)  ("  S2  +  ^4) 

Ah  =  mi  (~  ®i  “  ^2  “  S3  +  ^4  "  S5  +  ^7  +  ^9)  +  (l  +  (-  S6  +  S8  +  s10) 

Acceleration  terms  due  to  precession  of  the  earth’s  equator  were  developed  but 
found  to  be  well  below  required  precision. 

Lunar  and  solar  tide  coupling  effects 

The  principal  component  of  orbital  acceleration  due  to  tidal  forces  is 
undoubtedly  in  the  transverse  direction.  We  will  assume,  without  validation,  that  any 
radial  and  normal  components  of  such  acceleration  produce  only  negligible  effects. 
Then  the  moon,  in  its  orbit  around  the  earth,  would  have  an  additional  term  of  the 
form 


he  X  rG 

=  cw 


The  semimajor  axis  of  the  moon  is  a  good  approximation  to  rG,  so  that  we  can 
use 


.  ~  „  ^8  X  ?G 

ArG  “  c  0.U0256  hG 


(21) 


This  is  how  we  felt  when  this  force  was  First  introduced.  Now  we  are  no  longer 
sure  that  putting  \  for  rG  is  a  satisfactory  approximation,  especially  when  (21)  is 
employed  later  to  obtain  partial  derivatives.  No  changes  in  the  model  were  made, 
though. 


vector 


In  (21),  as  in  many  other  parts  of  this  paper,  the  angular  momentum 
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h  =  r  X  r. 


For  the  reader  familiar  with  the  treatment  of  tidal  effects  in  general  meory,  the 
following  remarks  should  be  made.  Taking  Gauss’  equation  for  the  time  derivative  of 
the  semimajor  axis  and  considering  (21),  we  have 


a 


C 


Double  integration  gives,  in  mean  anomaly  or  longitude,  the  term 

Sl  =  ^(t-t0)2 


Back  to  the  problem  at  hand.  Using  (21),  we  need  to  determine  the 
corresponding  heliocentric  accelerations.  The  barycenter  of  the  earth-moon  system  is 

.  _  =  m3F3+mnFn 

rB  m3  +  m,  j 


The  same  holds  for  the  accelerations  and,  in  particular,  for  a  small  increment: 


(m3+mn)ArB  =  m3Ar3  +  mj  j  AFj  j 


Since  the  perturbations  of  tidal  forces  on  this  barycenter  are  insensibly 
small,  it  is  quite  safe  to  put 


afb  =0 


There  follows 
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(22) 


There  is  no  contribution  to  the  motion  of  the  moon.  Indirect  effects  are 
automatically  applied  through  the  equations  of  motion.  In  this  case,  replacing  r3  by 
a3  is  probably  an  adequate  approximation. 

The  Oblate  Sun 

It  will  be  assumed  that  only  the  second  zonal  harmonic  is  needed. 
Letting  z$  be  a  z-coordinate  in  the  equatorial  frame  of  the  sun,  the  disturbing 
potential  is 


(26) 


The  meaning  of  the  various  symbols  is  believed  to  be  clear.  Since  only  J2  is 
considered,  the  subscript  2  was  left  off.  We  used  the  number 


Rq  =  0.0046524  a.u. 


J0  will  be  solved  for,  and  an  initial  estimate  of  zero  will  serve  nicely. 

Before  computing  the  accelerations  due  to  (26),  we  only  have  to 
relate  zs  to  the  available  coordinates.  The  procedure  is  quite  analogous  to  that 
employed  in  reducing  the  lunar  equatorial  elements.  The  transformation  matrix  here 
will  be  called  N,  but  it  has  the  form  of  M  in  the  earlier  section.  Since  only  z$  is 
needed,  only  the  third  row  of  N  will  be  required,  labeled  N3.  Hence 


sin  £2  sin  I 
cos  S2  sin  I 
cos  il  sin  I  sin  e  +  cos  I  cos  e  > 


N3  =  I  -  cos  S2  sin  I  cos  e  -  cos  I  sin  e 


where 

S2  is  the  longitude  of  the  ascending  node  of  the  solar  equator  on  the  ecliptic 
and 
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I  is  the  inclination  of  this  equator  to  the  ecliptic. 


In  the  Explanatory  Supplement  on  page  307  we  find,  for  our  frame 

of  1950: 


SI  =  75°4' 
I  =  7° 15' 


There  follows  that 


/  0.12194  \ 
N3  =  [-0.42454  ] 
\  0.89716/ 


with  much  more  precision  than  required.  We  can  now  put 


zs  =  N3  •  r 


into  (26)  and  start  differentiating.  The  resulting  acceleration  terms  are 

k-sfo-rH 


3 

Ar  =  -  - - 


r  +  N3  •  r 


Note  that  p  =  k2  because  of  our  units. 


(27) 


Above  equation  is  valid  for  all  planets.  However,  since  only  the 
observations  of  Mercury  and  Venus  are  to  be  used  to  determine  Je ,  (27)  was  added 
only  to  the  equations  of  motion  of  these  two  planets.  This  is  more  than  adequate, 
as  may  be  seen  later. 
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4. 


INTERPOLATION  AND  LIGHT  TIME 


The  coordinates  and  numerically  integrated  partials  will  have  to  be 
interpolated  to  the  observation  time.  We  are  using  the  Lagrangian  interpolation 
polynomial  of  degree 


m  +  4  m  even 
m  +  3  m  odd 


This  m  is  meant  to  be  the  order  of  the  running  routine.  The  reasons  for  the 
relatively  high  degree  are  as  follows.  Since  the  accelerations  are  integrated  with  a 
polynomical  of  degree  m,  the  corresponding  coordinates  require  m  +  2  if  to  be 
interpolated  without  loss  of  accuracy.  But  since  our  formulation  accepts  only  even 
b,  and  because  of  considerations  involving  light  transmission  time  corrections,  we 
adopted  above  scheme. 

Until  now,  all  observation  times  were  light  front  arrival  times  at  the 
earth.  However,  the  coordinates  of  all  bodies  have  to  be  found  at  the  instant  the 
solar  illumination  reaches  them.  This  statement  also  implies  that  solar  observations 
are  the  only  ones  not  to  be  adjusted.  Hence,  for  moon  and  planets  one  first 
calculates 


t'  =  t- 


I  r j -  r3  | 


c 


where  c  is  the  speed  of  light. 

One  now  goes  back  to  the  interpolation  routine  and  finds  new  coordinates  Fj  at  t', 
but  always  retaining  r3  at  t.  The  process  is  continued  to  convergence,  that  is,  until 


It' 


n  +  1  tn 


I  <  T 
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For  us  r  =  10" 7  seemed  quite  satisfactory.  The  program  finally  interpolates  also  for 
velocities  for  the  last  t'. 

As  the  program  iterates  toward  a  least-squares  solution,  t'  will  also 
change  slightly  and  converge  to  a  final  value.  The  latter,  of  course,  is  not  known  in 
the  beginning.  We  investigated  this  problem  and  found  the  consequences  negligible, 
as  intuition  would  suggest. 

5.  RESIDUALS  AND  CORRECTIONS 

For  each  observation  we  make  the  computed  place 

*  =  tan"*  ^ 

x  "  x3  (29) 

z-  z3 

5  =  tan"1  =  ===== 

V(x-  x3)2  +  (y -  y3)2 


and  then  the  residuals 


Aa  (®OBS  ”  ®COMp)  COs5COMP 
■  ^OBS  ^COMP 


(30) 


Early  results  on  the  moon  indicated  a  distinct  break  in  the  declination 
residuals  at  about  JD  242  3859.0  (1924  March).  We  were  unable  to  find  a  probable 
cause  for  this  bias.  After  a  careful  determination  of  the  magnitude,  we  added  the 
correction 


55  =  +  0'.'57 


to  all  A5  j  j  before  above  date.  Note  that  this  adjustment,  together  with  the 
declination  bias  parameter  solved  for  of  about  0'.'33,  adds  a  large  0'.'90  to  those 
early  observations. 
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All  observations  except  those  of  Pluto  must  be  corrected  for  the 
effects  of  annual  aberration.  Hence,  we  add  to  their  residuals 

1,.  .  •  x 

6  a  cos  5  =  — (x  3  sin  a  -  y3  cos  a) 

55  =  “(x3  cos  a  sin  6  +  y,  sin  a  sin  6  -  z3  cos  6) 
where  c  is  the  speed  of  light. 

In  keeping  with  standard  procedures  of  reducing  photographic 
observations,  Pluto’s  observations  had  not  been  corrected  for  the  eliptic  terms  of 
aberration.  We  make  this  adjustment  now  by  adding  to  the  residuals  of  Pluto 

SacosS  =  0"061  cos  a  -  0'.' 33 6  sin  a 
55  =  -  0'.'061  sin  a  sin  6  -  0'.'33 6  cos  a  sin  5  +0'.'027  cos  5 


As  mentioned  elsewhere  in  this  paper,  there  exists  a  pronounced  effect 
of  the  phase  of  Venus  on  both  a  and  6  .  Corrections  are  made  at  this  point.  To 
the  above  residuals  we  add 

5a  cos  6  =  -0"9159-  0"02401  X  lO"2?  +  0"22745  X10"V 
65  =  -0"6789  +  0!'81930  XKT  V0'.'46772  X10“V  +  0'.'93085  X  10“  V 

The  angle  <p  is  the  difference  of  orbital  longitudes  in  the  sense  Venus  minus  Earth. 
It  is  given  by 


=  74?37  +  0?616515  [t(JD)  -  242  0000.5]  0  <  <  360° 


In  a  final  step,  a  test  is  made  to  ascertain  that  the  residuals  are  very 
small  numbers.  They  could  conceivably  be  in  the  vicinity  of  -360°  or  +360°.  If 
so,  360°  is  added  or  subtracted  in  (29)  before  scaling  by  cos  5. 
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6. 


PARTIAL  DERIVATIVES 


Partials  for  Planetary  Orbits 

The  orbital  parameters  originally  chosen  for  the  major  planets  are  the 
nonsingular  elements  devised  by  Cohen  and  Hubbard  (1962).  At  the  time  we 
constructed  our  computer  program  this  seemed  like  a  sensible  choice  mostly  because 
various  subroutines,  involving  these  variables,  were  available  and  carefully  checked. 
The  obvious  drawbacks  of  this  approach  are  a  more  complicated  program  since  the 
nonsingular  elements  need  to  be  related  to  the  classical  ones  or  position  and 
velocity.  The  frequently  used  partials  by  Eckert  and  Brouwer  would  serve  as  well. 

We  will  not  repeat  anything  that  is  easily  found  in  the  paper  by 
Cohen  and  Hubbard.  In  order  to  simplify  both  formulation  and  programming,  their 
variables  were  renamed  as  follows: 


Ej  q,_ i j 
Ec  =  ev 


E,  =  e 
6  y 


i  =  1,2, 3, 4 


As  a  precomputation,  the  initial  conditions  of  the  planets  have  to  be  converted  to 
the  Ej  at  t0.  All  necessary  equations  are  given  in  the  Cohen  and  Hubbard  paper. 
Also,  their  equation  (9)  has  the  form 


where  the  denominator  E^  +  E^  +  Ej  +  E^  has  been  absorbed  in  our  ajj.  We  compute 
and  store  away  for  future  use  the  first  six  ay.  The  third  row  is  not  needed. 

With  these  preparations,  we  are  ready  to  calculate  the  partials 


t 
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where  the  vector  D  contains  the  observables  a  and  8  and  E  our  six  nonsingular 
elements  Ej.  It  is  clear  that  the  values  of  these  elements  at  tQ  are  to  be  improved, 
but  we  will  skip  the  additional  subscript  zero.  Since  the  observed  a  and  8  of  any 
planet  are  affected  by  the  positions  of  planet  and  earth,  we  need  to  allow  for  this 
dependence  in  our  normal  equations.  Hence,  we  calculate  both  dtx  /  dE^  and 
0Dj/9E3  .  The  relation  between  D  and  E  is  easily  established  using  the  rectangular 
coordinates  as  an  intermediary,  that  is 


D,  =  D,[r f3(E,)] 


Therefore 


Hence,  (31)  can  be  written 


3Dj 

0D^  ^ 

1 

n 

drT  0Et 

d3  = 

9idA_ 

9E[ 

form. 

Because  of 

-,-T 

3r3 

drj 

3D; 

9A5. 

OJ 

mi 

to  i-J 

II 

drj  3Ej 

(30) 


(31) 


(32) 
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Equations  (30)  and  (32)  are  actually  employed,  and  evaluated  for  every  observation 
of  moon  and  planets.  In  case  i  =  10,  observations  of  the  sun,  only  formula  (32)  is 
calculated. 


The  first  member  on  the  right  of  (30)  and  (32)  can  be  obtained  from 
(28).  The  results  are 


(-Ay, Ax,  O' 
(Ax)2  +  (Ay)2  x 


_ 1 _ 

V(Ax)2  +  (Ay)2  [(Ax)2  +  (Ay)2  +  (Az)2] 


(-AzAx,-AzAy,(Ax)2  +  (Ay)2) 


where 


Ar  = 


V 


r3 


The  other  two  matrices,  namely  9?j  /  9E*  and  9r3  /  9E3  ,  are  given  by 
Cohen  and  Hubbard.  We  cannot  consider  the  second  matrix  a  special  case  of  the 
first  since  their  treatment  differs,  as  we  will  see  in  a  moment.  One  first  computes 


<P  =  n(t  -  t0)  + 


a2n 


Next  comes  p,  Cx,  and  CY,  given  by  Cohen  and  Hubbard  in  equations  (49),  (50), 
and  (51).  Then  we  need  their  quantities  X  and  Y  which  are  best  obtained  from 


X  -  at  j  X  +  3i2^  +  al3  z 
Y  =  a?  l  x  *  a22^  *  a23  z 
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Wherever  the  time  appears  in  <p,  pj  Cx,  and  CY>  it  is  to  be  taken  as  the  reduced 
time  t'.  When  dealing  with  the  partials  for  the  elements  of  the  earth,  recall  that  t' 
is  identical  to  t.  Finally,  the  partials  9r/9ET  are  given  by  the  six  matrices  in  the 
Cohen  and  Hubbard  equation  (48). 

The  above  equations  seem  to  suggest  that  all  elements  or  coordinates 
to  be  evaluated  at  the  running  time  t  or  t',  as  opposed  to  initial  values  at  tQ,  can 
be  taken  from  the  numerical  integration  results.  This  must  not  be  done.  The 
difficulties  encountered  when  evaluating  Keplerian  partials  using  ingredients  from  the 
actual,  perturbed  orbit  were  discussed  in  detail  in  an  earlier  paper  (Cohen,  Hubbard, 
and  Oesterwinter,  1967).  Hence,  we  proceed  as  we  did  in  the  1967  study.  For  the 
purpose  of  calculating  the  partials,  a  Keplerian  orbit  was  carried  for  each  planet 
which  coincided  with  the  actual  orbit  at  epoch.  This  is  simple  enough  since  only 
the  mean  anomaly  needs  to  be  updated  for  each  observation  time.  Subsequent 
conversion  to  coordinates  is  done  by  a  subroutine  already  available. 

Partials  for  the  Lunar  Orbit 

A  fair  amount  of  time  was  spent  in  experimenting  with  various  sets 
of  partial  derivatives  for  correcting  the  lunar  orbit.  It  is  clear  that  Keplerian  partials 
are  inadequate  because  of  the  very  large  perturbations  in  the  lunar  orbit.  We  did, 
however,  try  to  use  partials  that  took  account  of  the  secular  perturbations  in  1,  g, 
and  h.  This  approach  failed  miserably,  evidently  because  of  the  large  periodic 
perturbations.  We  then  appealed  to  lunar  theory  and  used  partials  considerably  more 
sophisticated.  Again,  no  reasonable  solution  was  obtained.  Since  such  analytical 
partials  are  used  successfully  elsewhere,  perhaps  we  carried  an  insufficient  number  of 
terms.  Op.  the  other  hand,  it  is  probable  that  the  addition  of  clock  corrections  and 
tidal  coefficients  made  for  a  fairly  ill-conditioned  normal  matrix.  In  that  case, 
partials  of  a  given  accuracy  might  fail  which  would  be  perfectly  good  when  solving 
for  orbital  parameters  only. 

By  this  time  we  had  decided  to  pay  the  price  for  partials  obtained 
by  numerical  integration  We  have  used  this  device  for  many  years  in  application 
where  analytical  theories  do  not  exist  or  where  great  precision  is  required.  The 
latter  is  the  case  in  our  work  in  satellite  geodesy  and  various  other  studies  involving 
artificial  satellites.  Many  significant  figures  are  lost  in  inverting  the  normal  matrix 
which,  obviously,  requires  that  the  partials  have  several  more  figures  to  start  with. 

Since  we  are  dealing  with  additional  differential  equations,  usually 
lengthy,  that  are  to  be  integrated  numerically,  this  approach  is  always  costly.  In  the 
case  at  hand,  we  have  30  differential  equations  for  the  motion  of  the  10  bodies 


/ 
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involved.  It  will  be  seen  that  the  partials  for  the  lunar  orbit  corrections  require  18 
more  equations.  Hence,  computer  time  is  increased  by  50%  or  so. 


The  differential  equations  in  question  are  usually  called  variational  or 
perturbation  equations.  The  procedure  depends  on  the  fact  that 


The  derivatives  exist  and  all  functions  involved  are  continuous.  P  is  any  parameter  r 
depends  on.  The  right  hand  side  of  (33)  can  always  be  obtained  since  the  equations 
of  motion  are  known.  Cowell  integration  of  this  expression  then  yields  the  desired 
partials  3r/ 3P .  In  our  problem,  we  need  the  matrix  3F/3E1,  where  E  contains 
any  six  elements  describing  the  lunar  orbit.  Hence  the  18  different  partials 
mentioned  earlier. 

We  decided,  somewhat  arbitrarily,  to  differentially  correct  the 
geocentric  orbit  of  the  moon.  It  is  possible  that  a  set  of  heliocentric  parameters 
would  work  equally  well,  especially  when  rigorous  partial  derivatives  are  used.  So  we 
need  the  geocentric  acceleration  of  the  moon  which  is  obtained  from  the 
heliocentric  equations  of  motion  and 


rG 


The  result  is  short  enough  to  be  recorded  comfortably: 


-  k2(m3  +  m,  j ) -j  -  k2  -  k2  I  m  J—1'  J3-jj  -  " "-?U'jy) 

rG  \rii  r7  j;\  v  rj  ~  r3  i3  irrrn  >7 

2  k2(m3  +mi  i ^20 R1  r^*  [(5  if"  ~  ”  2zGP] 


(34) 


This  equation  is  seen  to  be  incomplete.  It  does  not  contain  relativistic 
or  tidal  terms  nor  spherical  harmonics  except  the  second  zonal  harmonic  of  the 
earth.  We  made  a  calculation  to  determine  the  contribution  of  relativistic  effects  on 
certain  partial  derivatives  for  Mercury  and  found  figures  like  1.0  X  10" 7  in  units  of 


the  leading,  term.  Hence,  we  concluded  that  it  would  be  safe  to  ignore  relativity  in 
partials  for  the  moon.  The  other  neglected  terms  were  not  checked  carefully.  They 
were  assumed  to  be  negligible  because  of  their  relative  size  in  the  equations  of 
motion.  In  any  event,  with  partial  derivatives  based  on  (34)  the  least-squares 
solutions  converge  rapidly,  making  any  gross  errors  in  this  area  unlikely. 

The  effects  of  small  changes  in  the  geocentric  orbit  of  the  moon  on 
the  orbit  of  the  earth-moon  barycenter  around  the  sun  are  minute.  We  can  safely 
assume  that  the  heliocentric  coordinates  are  invariant  when  correcting  the  geocentric 
orbit  of  the  moon.  This  simplification  can  be  used  to  our  advantage  immediately.  It 
is  clear  that  we  can  replace  r3  and  Fj  t  in  (34)  by 

...  mn  _ 

r3  rB  m3  +  mt ,  rG 


r 


1 1 


m- 


rB  +  m,  +m 


1 1 


where  rB  is  the  barycenter  mentioned.  The  latter  will  be  treated  as  constant  when 
differentiating  (34).  Since  the  planetary  coordinates  can  also  be  treated  as  constants, 
the  only  variable  left  in  (34)  will 'be  rG. 

We  frequently  use  the  abbreviation 


where  P  is  any  parameter  r  depends  on.  In  this  particular  section  P  is  any  one  of 
the  moon’s  orbital  elements,  but  it  will  later  be  seen  to  be  advantageous  to  perform 
the  formal  differentiation  of  (34)  without  assigning  a  specific  meaning  to  P.  The 
final  result  of  the  operation  is 
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=  9jb  9£g 
?k  "  3Pk  "  3fG  3Pk 


,  1  /3  T 

k2(m3  +m.  1)73  (p-fG?G_13 

*G  YG 


m3+mil 


-  -k2(m  +  ni  ) 


C  R2  f~, 
V'201V3 


rw  v«ii3  5 

Z  rG 


l I  ■  ')V '  t  (7|  '  ) +  10|^ 2PpT. 


(35) 


k2m 


1 1 


,n3  +IT11  1 


9  Itl:  r  3  _ 

'  j?,  I  r  -  r  P  I?  -  r3l2  (rj "  r3)(rj "  r3} 

j-y!.3  L 


k2m„  9  ’  rrij 

2  — — h- 


m3  +mn  j=i  lf|-fn  |3 
j^3 


j^  pu(2  (fi”PnXVfn>  -}3 


What  cannot  be  seen  is  that  (35)  is  incomplete.  The  addition  will  be  treated 
where  it  logically  arises,  namely  under  tidal  effects.  The  unit  vector  P  has  nothing 
to  do  with  the  parameter  Pk.  As  elsewhere  in  this  paper,  P  =  (0,0, 1  )T  .  I3  is  the 
identity  matrix.  An  expression  such  as  rrT  is  shorthand  for  the  matrix 


48 


It  is  important  to  remember  that  all  coordinates  in  (35)  are  presented 
in  the  equatorial  frame  of  the  earth. 


Except  for  an  additional  application  to  be  discussed  in  the  next 
section,  equation  (35)  is  used  for  the  orbital  elements  of  the  moon.  Hence  Pk 
stands  for  such  elements,  and  we  found  it  more  convenient  to  deal  with  equatorial 
elements  rather  than  the  customary  ecliptic  ones.  What  remains  to  be  done  before 
(35)  can  be  integrated  numerically  is  to  find  the  initial  values  of  the  various  fk. 
This  is  easily  done  using  the  Keplerian  relations  between  elements  and  coordinates. 
They  are  given  here  in  a  form  we  like  to  use,  but  many  variants  can  be  found  in 
the  literature. 


3F  _  ] 

3a  a 


3r  cos  E  +  e  /  r  \  sin  E  • 

—  - - ; - r-r+(l  +  — ) F 

3e  1  -  e2  \  p/  n 


3r  * 
gf  =  N  X  r 


z  sin  h 
-  z  cos  h 
y  cos  h  -  x  sin  h 


3?  _  r 
31  n 


(36) 


3f 

9g 


h  X  r  = 


y/T^W 


esin  Er — --  r 
\  a2n  , 


3r 

3h 
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(39) 


!£u  =  .  iEu  ?h. 
seJ  arj,  alj 


It  is  cle:  that  we  will  not  use  (39).  Geocentric  observations  of  the  moon  cannot 
be  used  directly  to  correct  the  heliocentric  orbit  of  the  earth.  There  is  an  indirect 
adjustment,  but  this  discussion  belongs  in  a  different  section.  Hence,  we  will 
evaluate  only  (38).  The  first  member  was  described  earlier,  and  the  second  is  the 
result  of  the  numerical  integration  just  discussed. 

Note  again  that  we  have  chosen  to  deal  with  equatorial  elements  for 
the  moon  in  contrast  to  ecliptic  ones  for  the  planets.  The  algorithm  for  the  moon 
was  developed  later,  and  it  was  noticed  that  going  to  ecliptic  elements  would 
introduce  unnecessary  complications. 

Partials  for  Lunar  Tide  Coupling 

An  earlier  version  of  our  program  contained  analytical  partial 
derivatives  for  the  tidal  coefficients  which  proved  to  be  inadquate.  There  is  probably 
a  simple  explanation.  It  was  seen  before  that  very  accurate  partials  were  needed  for 
the  lunar  orbit.  Since  these  elements  and  the  lunar  tidal  coefficient  are  solved  for 
simultaneously,  partial  derivatives  for  the  latter  must  be  of  comparable  quality. 
Hence,  we  turned  to  “rigorous”  partials,  again  obtained  by  numerical  integration. 
The  derivation  required  some  care,  and  we  felt  the  essential  steps  should  be 
recorded  here. 

What  we  need  is  3Dj  j  /  3C  which  can  immediately  be  expanded  into 


^ li  =  Hi  +!2ll  JUll 

ac  brj  ac  arjj  ac 

We  can  take  advantage  of  a  relation  used  previously,  namely 


(401 


aDjLi  =  SDjj  (41) 

afT  arf1 
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In  order  to  put  this  into  terms  of  geocentric  coordinates,  we  recall  - 


m 


11 


m. 


fb  F3  +  m3  +  mt j  fg  Fi 1  "  m3  +  mt j  fg  ^2) 

Since  we  assume  that  the  barycenter  is  not  affected,  differentiation  of  (42)  gives 


$ 

-mn 

3_!a 

m3  +  ml  1 

3C 

3F11 

“  m3 

3fg 

3C 

“  m3  +  n»i  i 

3C 

(43) 


Upon  substitution  of  (41)  and  (43)  into  (40),  our  required  partial  becomes 

3D1 1  _  3Dt  t  3rG 
3C  "  3fT  3C 


(44) 


The  factor  3Dj  j  I  is  quite  familiar  by  now.  The  second  factor  is  to  be 
obtained  by  numerical  integration.  Hence,  we  will  derive  the  corresponding 
variational  equation  next. 

Consider  the  acceleration  of  an  orbiting  body  in  the  presence  of  an 
additional  perturbation  such  as  tidal  forces.  It  is  clear  that 


r  =  r[r(r0,  F0,C,t),  r(r0,  r0,  C,t)*  C] 


Formally,  then, 


9F (F,  r;  C)  3F  (F0>  F0,  C,  t)  3j  (f  f  C) 

4-  -  -  4.  - — - - — - 

3fT  9C  9C 


(45) 


On  the  left  hand  side  we  indicate  that  a  partial  with  respect  to  C  is  needed  holding 
the  initial  conditions  and,  of  course,  the  time  fixed. 

In  the  following,  we  will  use  a  shorthand  notation  for  (45).  Also,  wc 
will  apply  above  relation  to  the  geocentric  orbit  of  the  moon.  Then  (45)  becomes 


fopW  3Fg  9rc  (  afG  3FG  3FG 

3C  a fl  ac  3fg  ac  ac 

A  large  part  of  the  matrix  3FG  /  brj.  is  already  available  and  given  by  (35).  Hence, 
we  now  have  to  supply  the  contribution  due  to  our  tidal  forces.  The  acceleration  is 
given  elsewhere  in  this  paper  as 


afg 


0.00256  hG 


The  result  of  differentiating  this  relation  is 


(47) 


3A£g 

K 


0.00256hr 


i 


FgFg 


^G  ^3  ^  +  ^3  ^  ^3  +  u2  (fGX  hG)(fc 

“G 


X,hr)T| 


T 


(48) 


Note  that 


o\  /o  -h,  h2\ 

0  =  [  h3  0  -h,  1 

1  /  \-h2  th,  0  J 
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The  next  matrix  in  (46),  namely  9rG  /9rJ  ,  stems  entirely  from  (47).  The  result  of 
this  operation  is 


0.00256hr 


f'o'j  "  V1)  jj  Ob  X  <0  (fG  X  I'o)  J 


(49) 


Thirdly,  the  last  vector  in  (46)  is  the  explicit  derivative  of  (47),  which  yields 


X  rG 


ac  0.00256 


(50) 


Adopting  the  1-  -  notation  introduced  earlier,  equation  (46)  now 


becomes 


/arG  9ArG\  9rG  ^  9rG 
‘3rT%fj/'%F^'  +  8C 


(51) 


These  are  the  three  new  differential  equations  to  be  integrated  numerically.  The 
vector  fc  is  generated  in  the  process  much  like  the  velocity.  The  other  factors  are: 


is  given  by  (35) 


9Arc 
3f G 


arG 

9C 


is  given  by  (48) 


is  given  by  (49) 


is  given  by  (50) 
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The  initial  values  of  fc  and  £c  are  both  0.  Results  of  the  numerical  integration  are 
substituted  into  (44),  and  this,  in  turn,  goes  into  the  normal  equations. 


It  seems  that  this  is  the  logical  place  to  update  equation  (35). 
Presently,  it  is  of  the  form 

r*  ■  ^  r„  <52> 

Because  of  the  introduction  of  tidal  terms  and  the  latter’s  dependence  on  the 
velocity,  (52)  must  be  augumented  to  read 


9rG  3ArG\_  3FG  . 

- + -  ?k  +  — 

3F5  /  a fj 


The  matrix  coefficients  are,  of  course,  the  same  as  in  (51).  The  initial  conditions  of 
i=k  are  the  same  as  before.  For  the  new  vector  £k  initial  values  are  obtained  from 
(37),  again  putting  t  =  tQ. 

Partials  for  Solar  Tide  Coupling 

As  we  saw  in  an  earlier  section,  the  result  of  a  constant  acceleration 
in  the  transverse  direction  is  a  quadratic  term  in  the  orbital  longitude  or,  to  be 
more  precise,  in  the  mean  anomaly.  Since  we  will  ignore  the  eccentricity  of  the 
earth  orbit,  this  perturbation  would  also  appear  in  the  argument  of  latitude,  that  is 

3  C0 

u  =  uQ  +n(t-  tQ)-  -  —  (t-  t0)2  (54) 


Figure  4  relates  u  to  the  coordinates  and  to  a  and  5 .  After  some  obvious 
intermediate  steps,  the  necessary  partials  are  found  to  be 
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Ecliptic 


FIGURE  4 

Relation  of  Argument  of  Latitude  n  to  Polar  and  Rectangular  Coordinates. 
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(55) 


9a  _  3  (t  -  t0  )2  cos2  ex3  Ax +  y3Ay 
9C0  2  a3  cos  e  Ax2  +  Ay2 


95  3  (t  -  tn)2  sin  e  cos  ex3(Ax2  +Ay2)-  (cos2  ex3Ay-  y3Ax)Az 

9C0  2  a3  cos  e  ^AjT2  +  Ay2"  (Ax2  +  Ay2  +  Az2) 


where  again 


These  relations  are  probably  more  compact  in  terms  of  a  and  6,  but  our  computer 
program,  at  this  stage,  has  more  convenient  access  to  the  rectangular  coordinates. 

If  any  perturbations  in  the  earth’s  orbit  can  be  obtained  from 
observations  of  the  sun,  they  must  also  be  visible  in  the  apparent  positions  of  some 
of  the  nearer  planets.  Hence,  (55)  is  executed  for  observations  of  the  sun  and  the 
major  planets.  The  computer  program  treats  these  bodies  identically  except  that 
r10  =  0  always. 

Partials  for  Solar  Oblateness 

In  a  later  part  of  this  paper  we  will  show  that  the  solution  for  the 
solar  oblateness  coefficient  is  attempted  simultaneously  with  the  planetary  orbital 
elements.  Hence,  approximate  analytical  partial  derivatives  suffice.  We  selected 
Brouwer’s  artificial  satellite  oblateness  theory  (1959)  to  supply  the  required  relations. 
As  stated  earlier,  we  will  consider  the  leading  zonal  harmonic  only  and  label  its 
coefficient  JQ. 


Required  for  our  purpose  are  the  partials  9D/9J0.  This  can  be 

expanded  into 


9D  9D  9r  9F  9EC 

_  s  s 

9J  9rT  drj  9EJ  9J 

s  s 


(56) 


The  first  matrix  on  the  right  is  well  known  by  now.  The  subscript  s  designates 
coordinates  and  elements  referred  to  the  plane  of  the  solar  equator. 
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It  is  the  last  vector  in  (56)  which  we  want  to  obtain  by  Brouwer’s 
theory.  For  this  purpose,  let 

E  =  osculating  elements  at  any  time 

E0  =  osculating  elements  at  tQ 

Eq  =  Brouwer’s  mean  elements  at  t0 

The  mean  elements  at  running  time  are  not  needed.  We  can  now  write 

E  =  e[e0(e"  j),t,j]  (57) 


Even  though  E(E0,  t,J)  is  not  explicitv  given  by  Brouwer,  this  is  no  stumbling 
block.  Using  (57),  we  can  compute  the  partial  derivative 


3E(Ey.  t,J)  _  3E(E0  t,J)  8E(En,tltJ)  9E0(Ey,J) 

3J  aj  aj  (5S) 

Inspection  will  show  that  the  partial  on  the  left  of  (58)  is  obtained 
by  differentiating  the  theory  since  the  constants  of  the  Brouwer  orbit,  E",  appear 
explicitely.  We,  however,  require  a  partial  in  which  the  osculating  elements  at  epoch 
are  held  fixed.  Hence,  we  solve  (58)  for  the  first  term  on  the  right  hand  side. 
Again  in  shorthand,  this  means 


3E  =  _  8E_  3E,  (59) 

aj  aj  alj  aj 

Since  we  are  presently  working  in  the  plane  of  the  solar  equator, 

E  =  Es. 

The  first  member  on  the  right  of  (59)  is  quickly  obtained.  It  can  be 
shown  that  only  secular  terms  need  be  considered.  We  looked  at  the  coefficients  of 
most  periodic  terms  and  found  the  largest  about  10" 4  times  the  size  of  the  secular 
effects  over  50  yf,ars.  Hence,  all  periodic  terms  are  ignored  at  this  step.  Also, 
secular  terms  can  safely  be  restricted  to  first  order.  Copying  directly  form  Brouwer, 
we  have 
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A1  =  -|7'  T?n(l  -  302)  (t  -  tQ) 

Ag  =  -^72n  0  "  59 2)  (t_  lo)  (60) 

Ah  =  - 3y2 nO  (t  -  tQ) 

where 


i?  =  V  1  -  e2 


6  =  cos  I 


The  difference  between  osculating  and  mean  elements  is  of  no  consequence  here. 
With  this  information  we  find 

3E  (Efl  =  _  3  nRo(t-^  /  o 
3J  ~4  a2  n4  /  0 

(  °0- 
\  1-50 

\  20 


The  second  term  on  the  right  hand  side  of  (59)  requires  some  care.  It  is  helpfui  to 
write  this  out  in  detail: 
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Since 


9E  3E0 


3ET  9J 


1  0  0  0  0  0\ 


0  10  0  0  0 


0  0  1  0  0  0 


0  0  10  0 


0  0  0  10 


9an 


9J 

9ep 

9J 


9J 


91o 

3J 

9go 

9J 


ooooi I  \  IT 


91  3  n 


(62) 


it  seems  that  this  secular  term  will  soon  exceed  all  the  periodic  functions  in  (62). 
Hence  we  proceed  to  reduce  (62)  to 


Second  thoughts  prompted  us  to  go  back  to  this  assumption  after  final 
results  were  obtained.  A  careful  numerical  evaluation  of  (62)  showed  that  for 
Mercury,  with  our  particular  set  of  initial  conditions,  the  term  retained  in  (63)  is 
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about  103  times  the  largest  one  neglected  at  the  end  of  the  data  span.  In  the 
case  of  Venus,  however,  the  term  retained  is  of  the  same  size  as  those  neglected. 
Fortunately  we  find  that  here  (3E/3Ej)(3E0/3J)  is  about  103  times  smaller  than 
3E(Eq)/3J.  For  both  planets,  then,  final  partials  are  good  to  about  three  significant 
figures.  Not  as  good  as  we  would  like,  but  probably  adequate  for  the  purpose. 

Pack  to  the  derivation.  Except  for  the  factor  t-  tQ,  which  appears  in  both 
terms  on  the  right  of  (59),  all  symbols  can  be  replaced  by  their  numerical  values. 
Since  the  semimajor  axis  has  only  short-periodic  perturbations,  3a<)  /3J  is  obtained 
from  that  very  expression  in  Brouwer’s  theory.  We  employed  the  following  values  in 
these  calculations: 


Mercury 

Venus 

Solar  Equator 

Inclination  to  ecliptic 

7°0' 

3°24' 

7°15' 

Node  on  ecliptic 

47°44' 

76° 14' 

75°4' 

Inclination  to  solar  equator  (derived) 

3°22' 

3°51' 

Semimajor  axis 

.38710 

.72333 

Eccentricity 

.20561 

.00682 

Mean  motion 

.071425 

.027962 

Radius 

.0046524 

Mean  anomaly 

324°9' 

272°  r 

Argument  of  perigee 

28°49' 

54°27' 

l 

i 


f 


Venus: 


3E_ 

3J 


io-5(t-t0) 


(65) 


These  represent  the  last  factor  on  the  right  of  (56). 

The  next  step  is  the  premultiplication  of  3E/3J  by 
9rs 


Since  the  first  three  rows  of  3E  /3J  are  zero,  the  first  three  columns  of  (66)  are 
not  needed.  With  the  aid  of  (36)  the  other  three  columns  become 

3r  _  1  . 

31  n 


3r  .  -  e  sin  E 

~  =  Fr  +  Gr  where  F  = - 

VT^ei 


r2  1 

G= - ~ZZT' 

a2.V  1  -  e2n 
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If  we  temporarily  let  the  elements  in  the  vectors  (64)  and  (65)  be  a,  0,  and  y, 
then  the  product  we  are  after  can  be  put  in  the  form 


ill 

3E*  9J 


10'5(t-  t0) 


0FF.+I? 


(67) 


The  subscript  S  was  added  to  3E  /  3J  at  this  stage  since  it  will  be  needed  at  once. 

According  to  (56),  we  now  premultiply  (67)  by  3r/3rJ.  If  the  coordinates 
referred  to  the  solar  and  terrestrial  equators  are  related  through 

fs  =  Nr 


then  the  matrix 


3r 


3fJ 


=  Nt 


Coinbining  the  last  three  equations,  we  obtain 


3r  3rs  3E  T 
- — -  — 5  =  N 

SrJ  3Ej  3J 


io‘s(t- 10) 


0FNr  + 


(68) 


Since  it  can  be  shown  that 


where  N3  is  the  last  row  of  N,  (68)  becomes 


3r  3FS  3E, 
3rJ  aif  3J 


10'5(t-t0) 
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The  remaining  steps  are  trivial.  We  recall  that  a,  0,  and  7  were  numbers. 
The  same  is  true  for  N3,  used  in  a  previous  section.  Also,  various  parts  of  F  and  G 
are  constant.  Hence,  we  substitute  numbers  wherever  possible.  It  becomes 
advantageous  to  introduce  a  few  new  symbols,  but  there  is  no  need  to  record  any 
relations  since  everything  is  easily  verified.  Finally  all  is  put  into  (56)  and  we 
obtain  3D.  3D- 

— -= - -  10“ 

where  3JQ  drj 


S(t-t0)[BiFi+  (q  +  D^  +  HjXrJ  (69) 


-.70590 
34.9470 
320.776 
/-. 20533 
(  .71504 

\- 1.5 11 07 


-.00235 
6.1823 
23.5894 
/-.02111 
I  .07351 
V-  .15534. 


Bj  =  bj  sin  E; 

D,=  d,r? 

The  subscripts  1  and  2  refer  to  Mercury  and  Venus.  Recall  that  all  but  the  first 
three  significant  figures  are  illusory. 

Partials  for  Clock  Corrections 

We  decided  to  solve  for  corrections  to  our  observation  clock  at  a  number  of 
equally  spaced  discrete  points.  Since  the  program  still  had  room  for  another  15  to 
20  unknowns,  we  arbitrarily  picked  the  15  dates  1912.5,  1916.5,  •••  1968.5  to  be 
thus  improved.  As  long  as  ET  was  used,  we  had  to  assume  at  least  two  points 
known.  Hence,  the  values  of  ET  at  1912.5  and  1968.5  were  held  fixed.  After 
switching  at  AT,  we  solve  for  corrections  from  1912.5  to  1952.5.  Since  AT  is  well 
known  from  laboratory  determination  from  1955  on,  we  do  not  solve  for 
corrections  to  the  points  1956.5  to  1968.5. 

The  principle  of  finding  our  corrections  is  a  simple  one.  Let 


o 

T  =  T°  +  e 
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T  =  TO  +  e 


where  T  is  the  corrected  time,  T°  the  initial  estimate,  and  e  the  increment  to  be 
solved  for.  Since  e  is  quite  small,  the  right  ascension  of  any  one  of  our  objects  can 
be  written  as 


a(T)  =  a(T°)  +  ae 


(70) 


The  same  holds  for  8.  Assume  further  that  the  e  for  any  observation  time  t  is 
given  by  the  four-point  interpolation 


e 


e_,  v p(ej-  e_j)  +  -p(p-  l)(e_2  -  e_,  -  e,  +e2) 


(71) 


where 


t-t-i 

1461 


and  1461  is  the  number  of  days  in  four  years.  The  subscripts  do  not  follow 
adopted  notation,  but  they  provide  a  more  symmetric  form.  It  is  almost  obvious 
that  t_  2  <  t_j  <  t  <  tj  <  t2-  The  t  and  e  with  subscripts  belong  to  grid  points 
for  which  corrections  are  computed. 


The  partial  required  is  da/de j.  This  can  be  written  as 


9a  9a  9e_ 
9e{  9e  9£j 


(72) 


The  first  factor  on  the  right  is  found  with  the  aid  of  formula  (70),  that  is 


9a 

9e 


=  a 
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From  the  relations  between  coordinates  and  a  and  5  we  obtain 


9e 

<JF7  =  0 


1-p 


de 

d£l 


=  P 


(75) 


To  summarize,  the  required  partials  are 


9D  _  _L_3e_ 

3ei  °9ei 


with  D  given  by  (73)  and  3e/0ej  by  (74)  or  (75). 

After  a  solution  for  the  e’s  if  made,  values  are  computed  for  the  years 
between  grid  points.  For  the  years  1913.5,  1914.5,  and  1915.5,  this  is  again  done 
by  linear  interpolation.  Four-point  interpolation  is  used  from  1917.5  on.  The 
complete  list  of  e’s  then  replaces  the  AT  -  table  in  our  Data  Conversion  Program, 
discussed  in  the  section  on  Data  Preparation. 

Partials  for  Lunar  Declination  Bias 

Once  lunar  residuals  came  down  to  a  reasonable  level,  it  became  apparent  that 
there  existed  a  severe  bias  in  declination.  The  same  has  been  seen  by  other 
investigators  of  lunar  observations,  and  the  cause  is  still  under  investigation.  We 
elected  to  “remove”  the  problem  by  solving  for  a  bias  parameter  which,  in  turn,  is 
applied  to  the  declination  residuals  of  the  next  iteration.  If  we  let  the  observed 
declination  be 


5  =  6°  +  A6 


where  A5  is  the  parameter  we  want,  then 


With  this  definition,  A5  is  added  to  the  .  ornputed  declination. 
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NORMAL  EQUATIONS  AND  THEIR  SOLUTIONS 


Let  us  start  with  the  equations  of  condition 


“3J,  AE  =  AD  (76) 

0ET 


As  before,  D  contains  the  observables  a  and  5.  E  is  any  parameter  we  solve  for, 
not  just  an  orbital  element.  The  normal  equations  are  then 


^W-^AE  =  ^ 
9E  0ET  0E 


WAD 


which  we  frequently  abbreviate 


BAE  =  AN 


(77) 


The  weiglit  matrix  W  contains  only  22  distinctly  different  elements,  one  each  for 
the  a  and  5  of  the  ten  bodies  observed,  with  an  extra  set  for  the  moon.  Hence, 
only  22  numbers  need  be  stored.  The  weights  themselves  are  obtained  as 


wa 


and 


where  aa  and  og  are  the  formal  standard  deviations  of  each  planet’s  residuals.  As 
such  they  would  be  numbers  obtained  in  the  previous  interation.  Their  derivation 
will  be  discussed  later. 


Since  in  most  problems  of  this  nature  the  number  of  observations  is 
large,  the  dimensions  of  (76)  are  enormous.  Hence,  it  is  impractical  to  evaluate  (76) 
and  do  (77)  numerically.  We  therefore  evaluate  (77)  at  once.  We  also  found  it 
advantageous  to  take  one  observation  at  a  time  and  complete  various  calculations 
before  going  to  the  next  one.  Let’s  look  at  this  now. 

The  partials  for  the  Various  parameters  are  available.  They  are  of  the 
form  da  /  0E  and  05  /  0E.  Since  the  residuals  in  a  are  scaled  by  cos  5 ,  we  will 
have  to  do  the  same  for  the  corresponding  partial  derivatives.  For  the  sake  of 
brevity,  let 
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3a 

a  =  —  cos  5 
3E 


d 


35 

3E 


Then  the  normal  matrix  becomes 


B  = 


l Wa a i  +diWdd,  a,Waa2  +d,Wdd2  •••  ai  Waap  +  dt  Wddp 
0  a2Waa2 +d2Wdd2  •••a2Waap +  d2Wddp 


0  *“apWaap+  dpWddp 


and  the  right  hand  side  is 


AN  = 


/ajWaAa  +d1WdA5N 


a2WaAa  +d2WdA6 


hWia  +  d  W.A6 

\p  a  p  d 


Since  all  this  is  done  for  one  observation  at  a  time,  no  subscript  is  needed  to 
identify  th*  observation.  Hence,  the  subscripts  attached  to  a  and  d  refer  to  the 
parameter  involved,  and  p  is  the  total  number  of  parameters.  Since  B  is  symmetric, 
only  (p2  +  p)/2  distinct  elements  need  be  computed  whiie  the  others  are  supplied  by 
bjj  =  bjj  upon  completion.  It  is  important  to  recall  that  the  residual  Aa  already 
contains  cos  5,  as  discussed  in  an  earlier  section. 

While  the  above  numbers  are  available,  it  is  advantageous  to  perform  a 
few  other  calculations  at  this  point.  Hence,  for  each  observation  we  add  to  the 
sums 


2  W  (Aa)2  and  2  W„ 

a  a 
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After  processing  the  last  observation,  we  calculate 


S  Wa(Aa)2 
2  W. 


S/Na 


2  Wa(Aa)2 
% 


and  similar  expressions  for  the  declination.  Here  na  is  the  total  number  of  data 
points  in  a.  The  fir$t  gives  the  standard  deviation  of  a  single  a,  the  other  its 
signal-to-noise  ratio.  We  did  not  find  the  latter  to  be  useful.  Both  formulae  were 
evaluated  for  individual  planets  as  well  as  the  solar  system  as  a  whole.  Again  note 
that  a  contains  the  scale  factor  cos  8  which,  via  the  weights,  enters  B  and  AN.  • 

The  solution  of  the  normal  equations  (77)  is 


AE  =  B_1  AN 


(78) 


If  we  let  the  elements  of  the  inverse  matrix  B  1  be  br.}  the  formal  standard 
deviation  of  the  parameters  AE  in  (78)  are  given  by 


ai  “ 


(79) 


and  the  correlation  coefficients  between  parameters  by 


P 


ij 


o.n. 


(80) 


It  would  seem  that  AE  is  now  to  be  added  to  the  initial  estimates  of 
the  parameters.  We  found  quite  some  time  ago,  and  so  did  others,  that  this 
procedure  often  does  not  converge.  It  is  immaterial  whether  we  talk  about  our 
nonsingular  elements,  four  of  which  contain  the  semimajor  axis,  or  the  rectangular 
coordinates.  The  problem  arises  whenever  the  semimajor  axis  is  not  decoupled  from 
the  other  five  orbital  parameters. 
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It  is  not  difficult  to  see  the  reason  for  this  behavior.  Let  us  assume 

that  we  make  two  least-squares  fits  to  the  same  data  set,  solving  for  corrections 

first  to  the  classical  elements  A  =  (a,e,I,l,g,h)T  and  then  the  coordinates 

X  =  (x,y,z,x,y  ,z,)T.  The  notation  was  chosen  to  set  this  pheripheral  discussion  off 
from  the  rest  of  the  chapter.  Let  the  solution  of  the  former  be  A  A",  the  other  AX. 

Now  it  is  clear  that  the  crucial  part  of  the  first  solution  is  a  good  Aa.  Any  error 

in  a  will  produce  a  secular  effect  in  the  computed  places.  Since  this  is  not  so  for 
the  other  five  elements,  their  adjustment  is  less  critical.  We  have  to  make  the 

assumption  that  we  know  how  to  formulate  a  least-squares  algorithm  that  will  yield 

a  “good”  Aa  which  is  linear  in  the  residuals.  If  we  now  take  our  coordinate 

solution  and  determine  the  increment  in  a  associated  with  it,  we  find 


Aa  =  a(x+  AX)-  a(x) 


where,  of  course,  X  was  the  initial  estimate.  We  will  now  proceed  to  show  that  Aa 
can  be  quite  different  from  Aa,  and,  hence,  cannot  also  be  a  good  adjustment. 

A  Taylor  expansion  about  X  gives 


Aa  ~ 

9XT  2  axaxr 


AX  +••• 


Next  we  will  show  that  the  first  term  on  the  right  of  this  series  is  nothing  but  the 
good  Aa  discussed  above.  It  is  linear  in  the  residuals  because  AX  is.  We  can  start 
with  the  normal  equations  in  the  coordinates,  namely 


Since 


W  "TT"  AX  =  W  AD 

ax  axT  ax 


9D  ap  3A 

axT  3at  axT 
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this  becomes 


8AT  3DT 
3X  3A 


3AT  3DT 
3X  3A 


W  AD 


Eliminating  the  leftmost  matrix  on  both  sides,  there  remain  the  normal  equations  in 
the  elements.  Hence,  we  must  have 


3A  _  _ 

-=f  AX  =  AA 
3X 


or 


3XT 

Q.E.D.  If  we  can  also  show  that  the  second  and  higher  order  terms  in  above  Taylor 
series  can  be  significant,  then  Aa  will  be  nonlinear  in  the  residuals  and  poor 
whenever  this  occurs. 

As  an  example,  we  picked  an  intermediate  solution  for  Mercury.  The 
clement  solution  calls  for 


Aa  =  -.700  X  10" 8  a.u. 


Solving  the  same  normal  equations  for  coordinates  and  substituting  these  into  above 
Taylor  expansion,  we  obtain 


Aa  =  (-.700  +  .175)  X  10" 8  a.u. 


The  second-order  term,  namely  0.175,  is  relatively  larro,  but  not  large  enough  to 
prevent  convergence.  However,  in  our  example  Mercur  s  orbital  parameters  are  quite 
close  to  their  final  values.  In  fact,  the  A  a  of  .7  X  i0"8  a.u.  corresponds  to  only  8" 
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in  the  orbital  longitude  after  20,000  days,  some  227  revolutions  of  Mercury.  Since 
the  ratio  of  the  second  to  the  first-order  term  in  above  Taylor  series  is  proportional 
to  AX,  it  is  clear  that  this  ratio  can  exceed  unity  if  AX  is  large  enough.  In  our 

example,  if  AX  had  been  greater  than  .700/.  175  =  4  times  the  actual  AX,  Aa 

would  have  the  wrong  sign.  In  other  words,  if  we  had  an  initial  estimate  of 
Mercury’s  orbit,  at  our  particular  epoch,  with  an  error  of  about  3  X  10" 8  a.u.  in 
the  semimajor  axis,  any  attempt  to  differentially  correct  the  rectangular  coordinates 
would  lead  to  divergence. 

It  can  be  argued  that  choosing  the  mean  motion  as  one  of  the 

solution  parameters  would  be  superior  to  improving  the  semimajor  axis.  Indeed,  the 

right  ascensions  of  a  planet  are  linear  in  n  except  for  periodic  terms.  But  it  can  be 
shown  quickly  that  there  is  very  little  difference  between  solving  for  a  and  n. 
Suppose  a  solution  in  a  has  been  made.  Then  the  corresponding  correction  for  n  is 
obtained  from 


An  =  n(a  +  Aa)  -  n(a) 


Again,  by  Taylor’s  theorem, 


dn  1  02n  , 

An  =  —  Aa  + - (Aar  +  •  •  • 

da  ^  0a2 


For  our  example,  we  find 


An  =  (.194 +  .439  X  10“8)X  10“8 


Hence  the  second-order  term  is  down  by  eight  orders  of  magnitude,  a  very 
comfortable  ratio. 

Let  us  now  return  to  the  discussion  of  the  solution  process.  We  have 
seen  that  it  will  not  do  to  add  AE  to  the  initial  estimates  of  E.  Hence  we  will 
convert  AE  to  increments  in  the  classical  elements,  and  these  will  be  added  to  the 
old  values  to  get  the  new  ones.  We  will  go  to  the  classical  elements  via  rectangular 
coordinates,  but  this  intermediate  step  is  taken  only  for  convenience,  as  set  forth 
earlier.  We  now  look  at  this  sequence  of  operations. 


73 


/ 


Let  C  stand  for  position  and  velocity,  and  E  for  the  nonsingular 
elements.  Then  there  exists  a  relation 


C  =  C(E) 


from  which  we  obtain 


AC 


AE 


The  matrix  of  partials  will  be  labeled  M.  Hence 


AC  =  M  AE 


(81) 


permits  us  to  express  the  solution  in  the  nonsigular  elements  in  terms  of  equatorial 
coordinates.  Cohen  and  Hubbard  (1962)  furnish  the  elements  of  M,  as  explained 
under  Partial  Derivatives. 

In  order  to  get  the  variance-covariance  matrix,  rewrite  (78)  in  the 

form 


6DT 

AE  =  B" 1— -3-  W  AD 
3E 


A  similar  equation  exists  for  the  coordinates,  namely 


AC  =  B:1  —  W  AD 

c  ac 


These  expressions  are  related  through  (81),  and  we  find 


1  ,3DT  .  3DT 

B:1  -=r  =  MB'1  ~=T 
c  ac  3E 


(82) 
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Since  also 


D  =  D  [C(E)] 


there  follows,  after  differentiating  w.r.t.  E, 


3DT  _  _3C^  3DT 
9E  3E  3C 

Substitution  of  the  last  relation  into  (82)  gives  the  desired  result 


B;1  =  MB’1  Mt 


(83) 


This  equation  is  used  to  get  the  B" 1  -  matrix  for  the  coordinates.  Standard 
deviations  and  correlation  coefficients  are  computed  similar  to  (79)  and  (80). 

The  transformation  into  ecliptic  elements  can  almost  be  written  down 
from  inspection.  We  appeal  to  the  relation 


E 


E 


where  E  E  are  the  classical  Keplerian  elements  ir.  the  ecliptic,  and  CE  the 
corresponding  rectangular  coordinates.  Then 


a£e 


3Ee  3Ce 

acE  acT 


We  let  the  first  matrix  be  N,  the  second  eT.  Then 


AEe  =  NeT  AC 


(84) 


¥ 
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The  elements  of  the  inverse  of  N  are  given  by  (36)  and  (37).  Hence,  we  make  N-1 
by  using  these  relations  and  N  by  numerical  inversion.  Note  again  that  those 
calculations  are  made  at  t  =  tQ.  The  matrix  eT  is  also  available  and  defined  by 
(12).  It  is  now  easy  to  see  that 


Bg1  =  NeTB-leNT  (86) 

Again,  standard  deviations  and  correlation  coefficients  for  the  ecliptic  elements  are 
calculated. 


The  AEe  obtained  by  (84)  are  added  to  the  initial  estimates  of  these 
parameters.  This  step  completes  the  logical  loop  of  the  orbit  improvement  program- 
The  improved  Ee  are  now  converted  to  equatorial  coordinates  and  nonsingular 
elements.  This  is  the  procedure  for  the  elements  of  all  planets,  including  the  earth. 

Treatment  of  the  other  parameters  differs  somewhat.  Their 
contributions  to  the  matrices  M  and  NeT  are  identity  matrices  of  the  appropriate 
dimensions.  From  this  point  on,  all  but  the  lunar  elements  are  handled  almost  like 
the  planetary  orbital  parameters. 

In  case  of  the  moon,  we  obtain  ArG  and  ArG,  the  increments  in 
equatorial  geocentric  coordinates.  We  then  calculate 


rll,  FINAL  "  r3, NEW  +FG,  OLD  +_,  ,  _  AfO  ^86) 

II *  nij  | 

mn 

r3’  FINAE  '  f3>  NEW  ~m;+-m~  aFg  <87> 

and  similar  equations  for  the  velocities.  This  requires  some  explanations. 

The  coordinates  labelled  “OLD”  are  initial  estimates.  Those  with  the 
subscript  “NEW”  are  the  results  of  planetary  orbit  adjustments  discussed  just  above- 
One  of  the  functions  of  (86)  is  to  apply  the  same  corrections  to  the  heliocentric 
coordinates  of  the  moon  that  were  determined  for  the  earth  from  observations  on 
sun  and  planets.  In  other  words,  any  such  adjustments  are  made  to  the  earth-moon 
baryccnter.  This  step  improves  the  speed  of  convergence  considerably.  Moreover, 
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equations  (86)  and  (87)  also  show  that  the  corrections  found  for  the  geocentric 
orbit  of  the  moon  are  distributed  to  the  heliocentric  coordinates  of  earth  and  moon 
such  that  their  barycenter  remains  unchanged. 

It  is  not  necessary  to  solve  for  all  parameters  that  enter  the  normal 
equations.  The  program  permits  suppression  of  any  number  of  unwanted  parameters 
in  the  solution.  Suppose  the  kth  unknown  is  not  needed.  The  program  sets  the  kth 
row  and  column  of  B,  see  equatioion  (77),  equal  to  zero  except  for  a  1  in  the 
main  diagonal  position  bkk.  Also  the  kth  element  of  AN  is  put  to  zero.  The 
solution  then  proceeds  in  a  normal  fashion.  A  secondary  solution  program  allows  us 
to  try  as  many  combinations  of  suppressed  parameters  as  we  wish  to  test  with  a 
given  set  of  normal  equations. 


VII.  CHECKOUT  AND  EXPERIMENTS 


The  original  formulation  of  our  computer  program  was  checked  with  care, 
and  practically  every  phase  was  verified  by  independent  hand  calculations.  Later 
additions  and  modifications  were  perhaps  not  all  examined  with  identical 
thoroughness,  but  none  was  employed  without  some  kind  of  checkout.  Only  the 
important  or  interesting  results  of  checks  and  experiments  will  be  reported  in  the 
following  pages. 

1.  INTEGRATION  ROUTINE 

Although  known  to  be  of  limited  value,  we  numerically  integrated  a 
Keplerian  orbit  approximating  that  of  the  earth.  The  stepsize  was  determined  such 
that  the  orbit  would  close  after  460  steps,  simply  for  the  sake  of  convenience. 
Maximum  error  in  any  coordinate  was  5  units  in  the  last  place  of  the  14  digit 
word.  Velocity  components  did  even  better.  All  supporting  hand  calculations  were 
done  to  20  significant  figures. 

2.  TRUNCATION  ERRORS 

Initial  tests  were  made  when  our  plans  did  not  include  integration  of 
the  lunar  orbit.  Most  test  runs  were  done  for  the  pair  Mercury  and  Jupiter. 
Mercury,  of  course,  has  the  shortest  period  and  Jupiter  the  largest  mass.  But  there 
are  additional  considerations.  We  have  seen  in  several  other  applications  of  numerical 
integration  that  the  highest  frequency  present  in  a  system  will  determine  the 
truncation  error,  even  if  this  frequency  is  only  a  perturbation  of  small  amplitude.  In 
other  words,  planets  such  as  Venus  and  Earth  would  probably  require  the  same  step 
size  as  Mercury  simply  because  the  latter’s  frequency  appears  in  the  higher 
differences.  Also,  Mercury’s,  eccentricity  of  0.21  is  likely  to  produce  a  noticablc 
effect. 


We  made  very  extensive  tests,  varying  all  possible  parameters.  It  would 
go  beyond  the  scope  of  this  paper  to  describe  the  results.  Be  that  as  it  may,  all 
this  is  of  academic  interest  anyway  since  we  introduced  the  moon.  It  is  clear  that 
the  lunar  period  will  now  control  truncation  error,  and,  consequently,  step  size. 

Before  we  turn  to  the  moon,  a  few  observations  of  general  nature, 
from  the  Mercury-Jupiter  experiment,  should  be  recorded.  The  position  of  Mercury 
in  its  orbit  at  the  epoch  of  integration  is  quite  important.  This,  obviously,  is  a 
consequence  of  the  relatively  large  eccentricity.  Any  Truncation  errors  thus  produced 
can  be  minimized  by  the  proper  selection  of  the  starting  routine  controls  a  and  b 

t 
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which  are  described  above  in  NUMERICAL  INTEGRATION  ROUTINE. 
Unfortunately,  a  and  b  would  then  at  least  be  a  function  of  initial  phase  of  a 
planet  and  its  eccentricity.  This  is  a  complication  we  do  not  wish  to  pay  for. 
Further  tests  showed  that  a  =  b  and  a  +  b  =  c  is  the  optimum  choice  for  the 
general  case. 


Also  seen  in  earlier  studies  was  the  effect  of  the  initial  relative 

position  of  the  principal  perturbing  body.  This  remark  as  well  as  some  made  above 

applies  only,  of  course,  when  the  step  size  approaches  its  allowable  upper  limit. 

With  an  integration  interval  several  times  as  small  any  such  effects  disappear. 

A  number  of  considerations  led  us  to  the  decision  to  integrate  all 

bodies  with  the  same  step  size.  It  does  not  seem  that  we  pay  much  extra  for 
computer  time,  but  we  do  have  a  simpler  program. 

Let  us  now  look  at  the  truncation  errors  in  the  presence  of  the 
moon.  It  is  clear  that  the  maximum  errors  occur  in  the  earth-moon  system.  We 
soon  found  that  the  heliocentric  coordinates  of  the  moon  show  larger  errors  than 
those  of  the  earth.  For  convenience,  we  put  results  into  the  geocentric  frame  and 
convert  truncation  errors  in  the  coordinates  into  angular  measure.  For  an  integration 
over  57  years  (20,800  days)  with  different  step  sizes,  the  errors  found  are  given  in 
Table  II. 


Since  we  are  using  a  step  size  of  O'?  4,  the  error  of  0"08  seems  to  be 
considerably  in  excess  of  our  design  accuracy.  Fortunately,  the  least-squares  fit  will 
remove  a  good  portion  of  this  signal.  As  long  as  the  same  program  with  the  same 
order  and  step  size  is  used  in  subsequent  runs,  the  effect  on  the  initial  conditions  is 
of  no  consequence.  Hence,  the  errors  actually  encountered  arc  the  quadratic  and 
higher  degree  parts  that  cannot  be  absorbed.  Our  tests  show  that  the  maximum 
truncation  error  left  is  only  0"013.  Figure  5  shows  this  to  occur  about  8,000  days 
after  epoch.  The  plot  depicts  angular  truncation  errors  after  fitting  a  lunar  orbit 
integrated  with  h  -  0^4  to  synthetic  observations  which  were  generated  with  a 
much  smaller  step  size. 

However,  when  the  program  is  exercised  with  the  lunar  tidal 
coefficient  as  a  free  parameter,  the  truncation  error  signal  is  practically  absorbed 
altogether.  This  may  change  the  tidal  coefficient  by  about  1%.  Since  its  standard 
deviation  is  more  like  20%,  the  error  is  insignificant. 
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TABLE  II 


Maximum  Truncation  Errors  of  Moon  After  57  Years 


Step  size 

in  a  cos  5 

in  5 

0?8 

550" 

180" 

0.5 

1.2 

0.9 

0.4 

0.08 

0.06 

* 
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COMPARISON  WITH  OTHER  INTEGRATIONS 


Since  our  integration  routine  evolved  from  other  in-house  Cowell 
algorithms,  it  was  carefully  compared  with  the  “best”  of  these  programs.  Naturally, 
any  such  checkout  is  concluded  only  after  all  discrepancies  are  removed  or  well 
below  the  threshold  of  significance.  If  the  same  computer  is  used,  such  differences 
can  often  be  pushed  down  to  the  last  significant  figure  or  two  because  of  many 
similarities  in  the  algorithms. 

Comparison  with  programs  developed  elsewhere  is  usually  a  more 
convincing  exercise.  However,  there  are  always  differences  in  computers  and  the 
mathematical  models  that  have  to  be  considered  or  circumvented.  We  made 
arrangements  with  O’Handley  and  Mulholland  at  JPL  for  such  tests.  At  that  time, 
their  programs  ran  on  the  IBM  7094  and  Univac  1108,  with  all  or  most  calculations 
in  double  precision.  Also,  JPL  uses  two  separate  programs  to  handle  the  major 
planets  and  the  moon.  In  the  former,  the  earth-moon  barycenter  is  integrated 
numerically  and  the  motion  of  the  two  bodies  about  this  point  is  supplied  through 
additional  input,  not  generated  internally.  The  lunar  program  contains  the  planets,  of 
course,  as  perturbing  bodies,  but  their  orbits  are  not  subject  to  correction.  There 
were  differences  in  the  relativistic  terms,  oblateness  perturbations,  and  other  small 
effects.  It  was  relatively  easy,  however,  to  adjust  our  program  to  simulate  the  two 
JPL  routines. 


Results  of  the  planetary  test  portion  are  quite  satisfactory. 
Experiments  were  run  over  12  years  with  step  sizes  from  O'?  25  to  1  ^  00.  However, 
the  larger  step  sizes  are  likely  to  reflect  the  differences  of  the  algorithms.  Since  this 
facet  is  of  no  interest  here,  we  will  restrict  our  comparison  to  h  =  O'?  25.  Table  III 
contains  the  principal  results.  The  differences  in  coordinates  have  been  converted  to 
angles  in  seconds  of  arc,  third  column,  as  seen  from  the  earth.  We  did  not 
investigate  the  relatively  large  coordinate  difference  of  the  earth.  This  number  may 
be  due  to  some  modeling  difference  we  failed  to  catch.  Even  if  correct,  it  would 
increase  some  of  the  angular  discrepancies  by  a  factor  of  ten  and  would  still  leave 
excellent  agreement  for  the  purpose  at  hand. 

Comparison  of  the  lunar  programs  did  not  fare  as  well.  In  one  run 
for  4,400  days  (12  years)  at  h  =  0^  25  the  final  position  difference  between  NWL 
and  JPL  was  2.2  X  10" 8  a.u.  This  corresponds  to  an  unacceptable  angular  error  of 
l!'8.  Material  for  additional  tests  did  not  become  available.  However,  we  feel 
confident  about  our  routine.  Since  we  have  only  one  program  in  which  the  moon  is 
but  one  of  n  planets,  the  various  tests  discussed  above  check  planetary  and  lunar 
orbits. 
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TABLE  III 


NWL  vs.  JPL  Coordinates  After  4,400  Days 


Planet 

|Ar| 

Angular  measure 

Mercury 

6.9  X  10“2  a.u. 

!'2X  10" 5 

Venus 

8.2 

.6  X  10"  5 

E-M  Barycenter 

24.9 

- 

Mars 

2.4 

.9  X  10'6 

Jupiter 

2.8 

.1  X  1G'6 

Saturn 

2.9 

.7  X  10“7 

Uranus 

6.3 

.7  X  10“  7 

Neptune 

10.8 

.8  X  10“  7 

Pluto 

4.2 

.2  X  10“  7 
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PHASE  EFFECTS  FOR  VENUS,  MERCURY,  AND  MOON 


Even  cur  earliest  solutions  indicated  that  there  was  a  problem  with 
the  Venus  residuals.  The  root  mean  squares  of  residuals  in  both  right  ascension  and 
declination  were  larger  than  those  of  other  planetary  transit  circle  observations. 
Closer  examination  quickly  revealed  a  signal  with  a  585  day  period  or  so,  clearly 
the  synodic  period  of  Venus.  This  suggested  that  the  phase  corrections  applied  by 
the  observers  did  not  adequately  reduce  observations  *o  the  center  of  the  planet. 
Figure  6  shows  a  residuals  over  a  bit  more  than  a  synodic  cycle  and  also  a 
quadratic  fitted  to  the  combined  material  of  several  cycles  spread  over  the  entire 
20,000  days.  The  polynomial  is  given  in  the  section  on  residuals.  It  is  no  surprise 
that  correction  of  residuals  in  all  subsequent  runs  considerably  improved  the  r.m.s. 
The  discontinuity  of  about  3"  occurs  at  inferior  conjunction,  as  one  would  expect. 
Residuals  in  a  were  also  examined  for  their  dependence  on  6,  but  no  relation  was 
found. 


There  is  also  a  distinct  phase  correction  in  declination,  although  it  is 
not  as  spectacular.  Nevertheless,  a  polynomial  of  third  degree  was  needed  to  describe 
this  effect  within  comparable  tolerances. 

The  residuals  of  Mercury  do  not  show  problems  of  this  magnitude  so 
that  we  did  not  make  a  search  for  phase  dependence.  However,  it  is  likely  that 
such  an  additional  phase  correction  can  be  found. 

For  the  moon,  all  a  residuals  for  16  consecutive  synodic  cycles  were 
examined,  and  nothing  of  any  consequence  to  this  study  was  found.  There  is  some 
indication  of  a  possible  discontinuity  near  new  moon  in  the  same  senrc  as  observed 
for  Venus.  If  real,  the  effect  would  be  small.  We  applied  no  corrections. 

5.  OBLATENESS  PERTURBATIONS 

The  effects  of  the  leading  zonal  harmonics  of  the  earth  on  the 
motion  of  the  moon  can  be  assessed  easily  using  one  of  the  existing  general 
theories.  One  finds  that  the  largest  secular  perturbation,  that  in  the  argument  of 
perigee,  amounts  to  about  700”  after  20,000  days.  The  other  two  secular  !erms  are 
about  half  that  size.  Coefficients  of  periodic  terms  can  be  found,  too,  but  they  are 
very  small. 


It  takes  a  more  powerful  theory  to  obtain  coupling  between  oblatencss 
and  solar  terms.  We  found  one  such  relatively  large  and  important  term  by 
coincidence,  investigating  a  different  problem.  It  is  shown  in  Figure  7  and  represents 
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FIGURE  6 


Venus  Residuals  in.  Right  Ascension  Uncorrected  for  Phase  Effects 


the  effect  of  suppressing  J2  of  the  earth  in  the  presence  of  solar  perturbations.  Its 
amplitude  is  about  7"  and  it  has  the  nodal  period  of  about  18  years. 

The  leading  zonal  harmonic  of  the  moon  produces  effects  which  are 
considerably  smaller.  The  secular  perturbation  in  g  is  only  about  1 1"  after  20,000 
days.  The  new  feature  is  that  here  the  precession  of  the  lunar  equator  must  be 
taken  into  account.  According  to  Cassini,  it  has  the  same  rate  as  the  node  of  the 
orbit.  Hence,  the  lunar  equatorial  plane  completes  one  precession  cycle  in  18.6 
years.  The  effect  on  the  lunar  a  residuals  is  mostly  secular,  something  like  6"  after 
20,000  days. 


Partial  derivatives  of  right  ascension  and  declination  with  respect  to 
the  solar  oblateness  coefficient  J,  were  tested  for  a  series  of  observations  on 
Mercury  and  Venus  against  finite  difference  ratios.  It  seems  that  the  partials  for 
Mercury  usually  agree  to  about  two  significant  figures.  This  is  by  no  means 
spectacular,  but  we  still  hope  that  they  serve  in  the  very  limited  role  they  have. 
Venus  partials  are  sometimes  10-15%  off,  but  they  enter  with  considerably  less 
weight. 


6.  LUNAR  TESSERAL  HARMONICS 

The  second  degree  terms  with  coefficients  c22  and  s22  were 
introduced  into  our  program  at  a  rather  late  stage.  Admittedly,  these  additions  were 
made  without  careful  analysis  of  the  size  of  the  resulting  perturbations  but  rather 
prompted  by  the  thought  that  the  earth  always  remains  near  the  first  meridian  of 
the  moon. 


The  first  test  run  immediately  indicated  a  problem  since  the  right 

ascension  residuals  increased  rapidly  with  time.  It  looked  like  a  growth  cubic  in 

time  which,  extrapolated  to  20,000  days,  would  amount  to  several  degrees.  However, 
any  growth  in  a  other  than  linear  must  be  erroneous  since  the  total  work  done  by 
such  an  acceleration  must  be  zero,  taken  over  a  long  enough  time.  There  could  be 
such  an  effect,  of  course,  if  the  earth  were  to  oscillate  about  the  first  meridian 

with  a  bias.  Hence,  differential  corrections  were  made,  with  the  aid  of  finite 

difference  ratios,  to  the  initial  orientation  of  the  moon  and  to  its  orbital  rate.  The 
first  is  given  by  the  argument  of  latitude  of  the  first  meridian  at  epoch,  the  second 
must  be  the  draconic  month.  We  found  that  our  residuals  called  for  an  increase  in 
said  argument  of  0°.2 9  and  a  decrease  in  the  period  by  0?40  X  10~4/day.  Since  the 
latter  amounts  to  only  0?80  after'  20,000  days,  both  numbers  sound  reasonable. 
However,  all  that  was  accomplished  was  a  change  from  a  signal  cubic  in  time  to  a 
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quadratic  one,  stili  leaving  a  residual  of  -150"  in  Aacos5  after  20,000  days.  It 
should  be  noted  that  the  real  problem  is  iarger  than  the  150"  indicate  since  much 
of  'he  signal  has  been  removed  by  the  simulataneous  corrections  to  the  lunar 
elements. 


Careful  checks  did  not  reveal  any  errors  in  either  formulation  or 
computer  program;  they  appear  to  be  well  hidden.  We  finally  decided  to  bypass  the 
tesseral  harmonics  altogether.  It  can  be  shown,  fortunately,  that  the  effects  of  this 
step  are  inconsequential.  The  neglected  accelerations  are  about  1.5  X  10” 9  times  the 
central  term.  Since  we  are  dealing  with  an  almost  constant  force,  the  missing 
accelerations  can  be  simulated  by  changing  the  combined  masses  of  earth  and  moon 
by  a  like  percentage.  Note  that  this  sum  of  masses  is  known  to  only  six  or  seven 
figures.  However,  since  we  do  not  adjust  the  masses,  the  effects  of  the  missing 
terms  must  be  absorbed  by  some  other  parameter.  It  cannot  be  the  mean  motion 
be-v.use  of  the  enforced  fit  to  observations.  That  leaves  the  semimajor  axis.  The 
relative  error  da/a  is  then  aboui  4.5  X  10" 9  which  is  much  smaller  than  the 
computed  standard  deviation. 

7.  PARTIALS  FOR  TIDAL  EFFECTS 

These  partials  were  tested  in  two  ways.  Both  lunar  and  solar  tide 

partial  derivatives  were  compared  with  results  from  finite  difference  ratios.  Since 
those  for  the  moon  are  obtained  by  numerical  integration,  they  are  expected  to  be 
quite  accurate.  Agreement  is  usually  found  to  three  or  four  places.  In  case  of  the 
sun,  partials  are  analytic  approximations.  Here  a  number  of  discrepancies  up  to  7% 
were  noted.  In  view  of  tne  marginal  value  of  this  particular  parameter,  no  other 
tests  were  made. 

The  lunar  partials  were  subjected  to  an  additional  and  more 

comprehensive  test.  Synthetic  observations  were  generated  without  tidal  terms.  Next 
an  orbit  was  fitted  to  these  data  by  improving  only  the  tidal  coefficient.  Its  initial 
input  value  was  -0.80  X  1 0“ 1 3 ,  and  this  number  was  driven  down  to  -0.76  X  10'17 
in  three  iterations.  In  the  first  iteration,  residuals  were  about  9"  at  1,000  days 

from  epoch,  which  came  down  to  0"007  in  the  third  iteration.  This  rate  of 

convergence  instills  some  confidence  in  the  procedure. 
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VIII.  CONSTANTS  AND  NUMERICAL  RESULTS 

In  order  to  keep  all  numerical  data  together  in  one  place,  we  are  listing  the 
input  constants  as  well  as  our  results  ir.  this  chapter. 

The  planetary  masses  were  changed  several  times  during  the  course  of  this 
study.  The  set  finally  adopted  and  given  in  Table  IV  is  identical  to  that  used  by 
O’Handley  in  his  Development  Ephermis  69,  but  the  decision  to  go  with  these 
particular  numbers  was  somewhat  arbitrary.  O’Handley’s  figures  for  the  terrestrial 
planets  benefit  from  electronic  observations. 

More  recent  mass  solutions  exist,  based  on  more  electronic  data,  but  they  are 
in  a  constant  state  of  flux  and  agreement  between  the  various  results  is  not  at  all 
satisfactory.  Also,  several  of  these  determinations  employ  procedures  that  we  have 
investigated  and  found  to  be  inadequate. 

Other  constants  used  are  as  follows: 


k  =  0.017  2020  9895 
c  =  299  792.5  km/sec 
1  a.u.  =  149  597  900  km 

C20  =  -.001  082  610 

C30  =  .000  0G2  539 

c20  =  -.000  207  5 


The  speed  of  light  and  the  a.u.  are  used  only  in  the  light  time  correction. 

In  an  earlier  part  of  this  paper  we  explained  why  the  program  fails  when 
attempting  to  solve  for  all  parameters  simultaneously.  What  happens  is  that  errors 
creep  into  the  lunar  element  corrections  large  enough  to  produce  residuals  of  50"  or 
so  after  55  years.  However,  it  is.  possible  to  iterate  successfully  to  a  least-squares 
solution  by  alternating  between  two  carefully  chosen  blocks  of  parameters.  Even 
then  one  must  proceed  with  care.  The  lunar  initial  conditions  arc  exceedingly 


89 


TABLE  IV 


Reciprocal  Masses 


Body 

Reciprocal  Mass 

Mercury 

5 

983 

000 

«• 

Venus 

408 

522 

Earth 

332 

945.561  925  441 

Mars 

3 

098 

700 

Jupiter 

1 

047.390  8 

Saturn 

3 

499.2 

Uranus 

22 

930 

Neptune 

19 

260 

Pluto 

1 

812 

000 

Moon 

27 

068 

03 

o 

<1 

OJ 

o 

o 

o 

1)  The  numbers  for  earth  and  moon  are  a  consequence  of  a  reciprocal  mass  of 
328  900.1  for  the  earth-moon  system  and  a  mass  ratio  of  81.301.  Fourteen 
significant  figures  were  retained  in  the  individual  masses  because  of  the  wordlength 
of  our  computer. 


sensitive  to  minute  changes  in  other  parameters.  Small  adjustments  of  planetary 
elements,  very  close  to  final  convergence,  can  easily  cause  signals  of  20"  in  the 
residuals  of  the  moon.  Hence,  it  should  be  no  surprise  that  the  last  three  iterations 
were  devoted  to  correcting  the  lunar  elements,  the  clock  corrections,  the  declination 
bias,  and  the  lunar  tidal  coefficient.  The  normal  equations  at  each  iteration  were 
subjected  to  a  variety  of  test  solutions.  Hence,  it  was  easy  to  ascertain  that  any 
additional  planetary  corrections  called  for  were  always  well  within  their  own 
standard  deviations. 

Table  XI  lists  the  number  of  observations  used,  weights,  and  the  standard 
deviations  (root  mean  square)  of  their  final  residuals.  Observations  in  right  ascension 
and  declination  are  now  counted  separately.  This  procedure  has  recently  been 
adopted  by  other  researchers  in  combining  optical  data  with  other  types  of 
observations,  such  as  radar  ranges.  The  weights  are  those  actually  employed  in  the 
final  run.  They  are  1  /  a2,  where  a  is  the  standard  deviation  of  the  residuals  in  the 
previous  iteration,  converted  to  radian  measure.  It  is  easily  verified  that  the  sigmas 
in  this  table  would  lead  to  very  similar  weights.  The  derivation  of  the  sigmas  was 
discussed  previously.  It  is  important  to  note,  however,  that  the  weight  of  all 

residuals  larger  than  4a  was  put  to  zero.  In  the  case  of  Pluto,  this  gate  is  probably 
too  narrow.  In  a  test  with  a  20a  cut  off,  the  sigmas  for  Pluto  went  up  to  i  "43 
and  1  "66,  in  a  and  5!  In  the  same  test,  the  a  for  all  other  planets  were  "01  or 
"02  larger  than  those  given  in  Table  XI.  This  indicates  that,  except  for  Pluto,  very 
few  residuals  fall  outside  tho  4a-band.  Hence,  4a  seems  to  be  a  comfortably 

conservative  cutoff. 

Among  the  principal  results  of  this  study  are  the  improved  orbital 
parameters.  All  adjustments  were  made  at  the  epoch  of  integration,  namely 

JD  242  0000.5.  We  will  give  coordinates,  Table  XII,  at  this  date  as  well  as  the 

tliree  following  future  dates: 


JD  242  0000.5  =  1913  August  21.0 
244  1200.5  “  1971  September  6.0 
244  1600.5  =  1972  October  10.0 


244  1000.5  =  1973  November  14.0 


Osculating  elements  will  also  be  listed,  but  only  lor  the  first  and  last  of  above 
times.  Formal  standard  deviations  are  given  for  coordinates  and  elements  at  epoch. 
We  wish  to  stress  the  fact  that  these  as  well  as  all  other  standard  deviations  quoted 

for  the  final  results  are  obtained  from  a  total  solution  for  all  parameters.  Such 

numbers  are  larger  than  the  sigmas  associated  with  a  partial  solution  and,  we  hope, 
more  realistic.  All  coordinates  are  tabulated  to  the  full  14  place  capacity  of  our 
computer  in  order  to  facilitate  using  our  numbers  as  initial  conditions  for  other 
runs.  The  elements,  on  the  other  hand,  are  truncated  in  keeping  with  their  formal 
standard  deviations.  Since  the  improvement  of  the  lunar  elements  was  done  in  the 
somewhat  unusual  equatorial  frame,  we  will  list,  at  epoch,  both  the  equatorial  and 
the  derived  ecliptic  elements.  See  Tables  XIII  and  XIV. 

In  using  any  of  our  figures,  the  following  must  be  kept  in  mind.  Data  given 

for  JD  242  0000.5  incorporate  the  adjustments  procuded  by  the  last  iteration.  All 

others  are  consistent  with  the  initial  conditions  to  the  last  run  and,  hence,  reflect 
the  next  to  last  improvement  cycle.  The  differences  are  of  academic  interest  only 
since  they  are  always  much  smaller  than  the  respective  sigmas. 

As  backup  to  the  published  tables  we  intend  to  save  our  tapes  which  contain 
coordinates  at  interval's  of  400  days. 

A  few  comments  should  be  made  about  the  standard  deviations  given  with 
our  coordinates  and  elements.  Table  XII  shows  that  both  y  and  y  for  Pluto  are 
considerably  weaker  than  the  other  coordinates.  It  is  interesting  to  speculate  by 
considering  the  geometry.  Since  Pluto  lies  approximately  along  the  y-axis  around 
epoch,  it  would  seem  likely  that  angular  observations  yield  less  resolution  in  this 
direction. 

It  is  not  difficult  to  see  why  some  sigmas  are  given  to  two  significant 
figures.  We  aimed  at  quoting  the  same  number  of  decimals  for  all  three  components 
ol  any  given  position  or  velocity  vector.  In  each  case,  the  smallest  sigma  determined 
the  number  of  decimals  retained. 

Up  to  four  significant  figures  are  given  in  certain  standard  deviations  in 
Table  XIV.  In  each  such  case  we  are  dealing  with  a  small  eccentricity  or  inclination. 
We  would  guess  that  the  nonsigular  elements  1  +  g  and  g  +  h,  respectively,  would 
have  sigmas  in  the  last  decimal  quoted.  Although  the  latter  values  can  be  calculated, 
we  truncated  our  table  simply  by  inspection. 

Table  XIV  also  shows  that  the  accuracy  of  the  semimajor  axes  deteriorates 
quickly  as  we  go  to  the  outer  planets.  From  Saturn  to  Pluto,  the  corresponding 
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sigmas  increase  by  three  orders  of  magnitude.  This  is  partly  explained  by  comparison 
of  their  periods  to  our  observation  span  of  55  years. 

In  the  case  of  Pluto,  it  made  no  sense  to  quote  1,  g,  and  h  to  the  same 
number  of  decimals.  The  in-plane  elements  1  and  g  also  suffer  from  the  relatively 
short  data  span.  The  sigma  for  h,  on  the  other  hand,  shows  that  the  node  is 
sharply  defined,  thanks  to  the  substantial  inclination. 

The  ecliptic  elements  for  the  moon  are  derived  quantities.  The  principal 
least-squares  solution  was  made  in  terms  of  equatorial  elements,  given  in  Table  XIII. 
Again,  it  is  seen  that  the  large  equatorial  inclination  leads  to  a  much  better  defined 
node  than  the  ecliptic  counterpart. 

We  did  not  think  it  necessary  to  list  all  correlation  coefficients  of  our 
73  X  73  matrix.  In  fact,  since  the  computer  prints  coefficients  for  the  nonsingular 
elements,  equatorial  coordinates,  and  ecliptic  elements,  there  is  a  total  of  about 
8,000  numbers.  We  decided  to  tabulate  only  the  correlation  coefficients  between  the 
six  elements  of  any  given  planet  and  all  others  exceeding  about  0.10.  Table  XV 
gives  the  former.  It  is  arranged  to  conserve  space,  but  it  must  not  be  misread  to 
contain  correlation  coefficients  across  planets. 

Table  XVI  gives  all  planetary  elements  which  show  correlation  coefficients 
larger  than  0.10  with  the  lunar  elements  and  the  clock  corrections.  They  are  seen 
to  be  the  elements  a  and  1  of  the  terrestrial  planets  and  a  for  Jupiter.  C  is  the 
lunar  tidal  coefficient.  Note  that  the  lunar  elements  I  and  h  are  missing  because  of 
very  small  coefficients.  The  1  for  Venus  was  left  in  order  to  complete  the  upper 
portion  of  the  table. 

Practically  all  remaining  planetary  coefficients  greater  than  about  0.10  are 
listed  in  Table  XVII.  About  half  a  dozen  were  ignored  so  that  we  would  not  have 
to  add  another  line  for  just  one  number.  The  largest  of  these  was  a  0.13  between 
the  a:s  of  Mars  and  Jupiter. 

Finally,  Table  XVIII  gives  the  complete  set  of  correlation  coefficients 
between  the  tidal  term  C,  the  lunar  elements,  and  the  clock  corrections.  Some  of 
the  large  numbers  are  not  unexpected. 

The  only  parameter  not  mentioned  at  all  so  far  is  the  declination  bias  AS. 
We  found  it  only  very  weakly  correlated  with  anything  else. 
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One  of  the  principal  results  of  this  study  is  the  extrapolation  of  the  atomic 
time  scale  back  to  1912.5.  The  final  figures  are  given  in  Table  XIX.  The  column 
labelled  AT0  is  taken  from  Brouwer’s  work  (1952)  through  1948.5.  From  1949.5 
on.  these  ^umbers  are  obtained  from  the  American  Ephemeris.  The  next  column 
contains  our  solution  “atomic  time  minus  universal  time”.  The  symbols  AT  and  UT 
were  chosen  purposely  different  from  the  abbreviations  A.l,  U.T.,  and  U.T.l,  which 
have  a  very  specific  meaning.  As  stated  elsewhere,  AT-UT  was  solved  for  at  intervals 
of  four  years,  and  intermediate  values  were  iupplied  by  four  point  interpolation. 
The  standard  deviations  listed  are  again  obtained  from  a  global  solution.  In  a  partial 
solution  for  clock  corrections  and  lunar  elements  only,  some  of  the  sigmas  were  only 
30%  of  the  tabulated  ones. 

The  final  column  in  Table  XVII,  AT-ET,  is  simply  (AT-UT)  -  ATQ.  It  is  a 
measure  of  the  disparity  between  the  atomic  and  ephemeris  time  scales.  Figure  8  is 
a  plot  of  these  numbers.  While  there  is  a  total  excursion  of  six  seconds  or  so,  the 
fine  structure  on  the  left  is  well  within  the  noise.  However,  some  of  the 
fluctuations  in  the  right  hand  half  of  the  figure  seem  reasonably  well  established. 

The  last  two  parameters  of  our  principal  solution  are  the  lunar  tidal 
coefficient  and  the  declination  bias.  For  better  visibility,  they  are  put  in  Table  V. 
The  product  KT2  would  be  added  to  the  computed  mean  longitude  of  the  moon, 
with  T  measured  in  centuries.  Over  our  data  span  of  55  years,  this  effect  amounts 
to  a  sizeable  6". 

The  bias  term  A8  of  -0'.'33  is  to  be  added  to.  the  computed  declinations  of 
the  moon  in  order  to  match  observations.  These  computed  places  are,  of  course, 
dynamically  consistent  with  all  other  ephemerides  of  our  solar  system.  Most  of  the 
observed  declinations  are  those  recently  uniformly  reduced  (Adams,  Klock,  and 
Scott,  1969)  and,  in  particular,  adjusted  using  Watts’  limb  corrections.  The  data 
before  1925,  not  treated  by  above  authors,  are  effectively  in  the  same  declination 
system  since  we  had  added  a  constant  bias  correction  of  0'.'57  to  these  early 
residuals.  Perhaps  the  center  of  Watts’  reference  sphere  is  0"33  south  of  the  center 
of  mass. 

Figures  9  to  28  are  plots  of  the  residuals  in  right  ascension  and  declination. 
In  viewing  these  it  should  be  recalled  that  all  points  greater  than  4o  have  been 
deweighted  and  do  not  contribute  to  the  solution.  The  Venus  residuals  still  seem 
too  noisy,  but  no  additional  systematic  signal  was  found.  From  Mars  to  Neptune, 
some  plots  exhibit  signals  which  appear  to  be  above  the  noise  level.  Clearly,  these 
trends  were  not  removable  by  improving  the  initial  conditions.  Whether  mass 
corrections  can  be  extracted  remains  to  be  seen. 


FIGURE  8 

Atomic  Minus  Ephemeris  Time,  Extrapolated  Back  to  1912.5 
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TABLE  V 


Tidal  Coefficient  and  Declination  Bias 


'i  C 

Note:  K  =  -  -  — 
2  a 


(-19"  ±  4")  /  century2 


-"33  ±."01 


/ 


Table  VI  contains  two  parameters  which  were  not  part  of  our  principal 
solution.  The  numbers  shown  here  stem  from  various  experimental  solutions,  but 
they  are  still  interesting. 

Kq  results  from  modeling  a  tidal  term  in  the  apparent  motion  of  the  sun  or, 
what  amounts  to  the  same,  in  the  equations  of  motion  of  the  earth.  Van  Flandern 
correctly  predicted  that  we  should  not  find  any  term  of  significance,  and  the  result 
bears  him  out.  is  seen  to  be  a  positive  quantity.  Since  we  are  not  aware  of  any 
forces  which  would  increase  the  earth’s  mean  motion,  we  omitted  this  term  from 
our  main  solution.  Since  the  numerical  value  is  3a,  it  has  some  statistical 
significance.  However,  the  total  effect  of  this  term  would  be  only  0'.'2  over  our  data 
span.  It  is  more  than  likely  that  there  are  a  number  of  unmodeled  forces  which 
could  give  rise  to  perturbations  of  this  size. 

The  other  quantity  in  this  table,  J0,  is  a  bit  of  a  disappointment.  We  had 
hoped  to  get  the  leading  zonal  harmonic  of  the  sun’s  oblate  field.  However,  since  it 
is  only  of  the  order  of  its  own  sigma,  it  has  really  no  significance.  Hence,  the 
wrong  sign  means  little.  If  JQ  were  actually  of  this  size,  it  would  produce  a  secular 
perturbation  in  Mercury’s  periheiion  of  about  2"  in  55  years,  and  about  1”  in  1  and 
h.  All  that  can  be  said  from  our  result  is  that  the  sun’s  J2  is  not  likely  to  be 
much  larger  than  1 .4  X  10~s . 


TABLE  VI 


Suppressed  Parameters 


K* 

(  +."6  ±  ."2)  /  century' 

4 

(—.17  ±  .14)  X  1CT4 
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IX.  SUMMARY  AND  RECOMMENDATIONS 


We  have  derived  and  presented  very  accurate  orbital  elements  for  the  moon 
and  the  major  planets.  It  is  our  belief  that  ephemerides  .based  on  our  initial 
conditions  are  a  measurable  improvement  over  numbers  presently  published  by  the 
international  astronomical  community.  In  fact,  our  solution  appears  to  be  almost  as 
refined  as  the  present  state  of  the  art  and  the  pocketbook  permit. 

Our  investigations  also  yielded  some  by-products,  and  two  of  these  have 
considerable  scientific  significance.  The  backward  extrapolation  of  the  atomic  time 
scale  extends  from  15  to  almost  60  years  the  span  over  which  a  clock  with 
constant  rate  is  now  available.  The  secular  deceleration  of  the  moon  has  been 
determined  over  the  same  interval.  Since  this  result  is  also  believed  to  be  a 

refinement  over  previous  numbers,  it  may  well  have  a  bearing  on  studies  of  the  past 
history  of  lunar  motion. 

In  the  introduction  of  this  report  we  have  explained  that  certain  deficiencies 
in  the  available  ephemerides  prompted  us  to  undertake  this  study.  Now  that  good 
results  are  in  hand,  our  initial  conditions,  or  those  stemming  from  another  modern 
global  solution,  should  be  made  the  basis  for  future  publications  of  the  positions  of 
sun,  moon,  and  planets.  We  believe  that  this  suggestion  merits  careful  consideration 
particularly  in  view  of  the  fact  that  the  International  Astronomical  Union  plans  to 
extend  the  present  set  of  tables  for  another  decade. 

At  the  moment,  we  feel  most  comfortable  with  our  own  results  and 

procedures.  As  explained  earlier,  the  simultaneity  of  solving  for  the  orbits  of  moon 

and  planets  is  a  crucial  feature  expecially  in  deriving  highly  accurate  lunar  elements. 

Also,  we  have  subjected  our  model  and  program  to  a  variety  of  tests,  and  they  now 

appear  to  be  free  of  detectable  errors.  Nevertheless,  we  wish  to  recognize  the 

somewhat  similar  efforts  of  our  friends  at  MIT  and  JPL.  Theirs  are  also  modern 

solutions  having  their  own  strong  points. 

The  introduction  of  ephemerides  computed  by  numerical  integration  and 
based  on  global  solutions  is  seen  as  a  first  step  only.  In  order  to  maintain  a 
product  that  is  to  be  useful  in  these  days  of  space  age  technology,  improvements 
must  be  made  continuously.  Different  and  more  precise  observational  data  will 

become  available  so  that  both  the  mathematical  model  and  the  computational  tools 
require  updating.  At  the  same  time  the  needs  of  interplanetary  navigation  will 
become  more  stringent.  Researchers  will  continue  to  find  that,  in  filling  the  practical 
requirements,  exciting  scientific  discoveries  are  made. 
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TABLE  VII 


USNO  Transit  Circle  Data  Punched  at  Dahigren 


Volume 

Instrument 

Time  Span 

Source  Codes 

XI 

6" 

« 

1911  - 1918 

D 

6 

XIII 

9" 

1913-1925 

D 

7 

XV  Part  V 

X-  g" 

1935-1944 

D 

8 

XVI  Part  I 

6" 

1925-1941 

D 

9 

XVI  Part  III 

6" 

1941  - 1948 

D 

0 

A-I 
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A-2 


TABLE  IX 


Identification  of  Pluto  Observations 


Source  Code 

Meaning  < 

Source 

D  PR 

Prediscovery  observations 

AJ. 

72,973  Table  1 

D  SB 

Sharaf-Budnikova  normal  places 

A.J. 

72,973  Table  II 

D  YM 

Yerkes-McDonald  observations 

A.J. 

72,973  Table  III 

D  LO 

Lowell  observations 

A.J. 

72,973  Table  IV 

D  YE 

Yerkes  observations 

A.J. 

72,973  Table  V 

D  HA 

Halliday  observations 

A.J. 

72,973  Table  VI 

DL  1 

Dahlgren  List  1 

Chernykh  and  Chernykh 

A-3 
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TABLE  X 


Adams  and  Cowell  Coefficients 


0 

+32,011,868,528,640,000 

+32,011,868,528,640,000 

1 

-16,005,934,264,320,000 

-32,011,868,528,640,000 

2 

-2,667,655,710,720,000 

+2,667,655,710,720,000 

3 

-1,333,827,855,360,000 

0 

4 

-844,757,641,728,000 

-133,382,785,536,000 

5 

-600,222,534,912,000 

-133,382,785,536,000 

6 

-  456,783,1 10,784,000 

-  116,974,585,728,000 

7 

-363,891,528,000,000 

-100,566,385,920,000 

8 

-299,520,219,398,400 

-86,707,632,211,200 

9 

-252,655,401,398,400 

-75,398,324,601,600 

10 

-217,227,737,563,200 

-66,193,573,118,400 

11 

-189,640,115,028,000 

-58,648,487,788,800 

12 

-167,636,336,098,320 

-52,401,453,198,480 

13 

-149,735,464,049,160 

-47,174,128,491,600 

14 

-134,928,496,929,540 

-42,755,108,505,900 

15 

-122,506,205,369,730 

-38,983,584,907,800 

16 

-111,956,703,448,001 

-35,736,323,456,205 

D  =  32,011,868,528,640,000 


All  numbers  are  to  be  divided  by  the  common  denominator  D.  Table  applicable  to 
all  orders  up  to  16. 
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TABLE  XI 


Number  and  Quality  of  Observations 


Body 

Number  of 
Observations 

1  Weight  X  10r11  in  Final  Standard  Deviation  a  in 
Aacos5  A5  Aac  os  5  AS 

System 

40  972 

■- 

*  • 

.65 

.72 

Mercury 

4  380 

.52995 

.59585 

.90 

.84 

Venus 

6  058 

.53712 

.50817 

.89 

.92 

Mars 

1  284 

1.42713 

1.39628 

.55 

.55 

Jupiter 

1  748 

2.06414 

1.56738 

.45 

.52 

Saturn 

1  750 

1.66827 

1.49200 

.51 

.53 

Uranus 

1  822 

3.21105 

1.76477 

.36 

.49 

Neptune 

1  796 

2.50643 

1.92599 

.4! 

.47 

Pluto 

1  110 

.54076 

.62207 

,892 

.83  2 

Sun 

14  148 

.74048 

.63428 

.76 

.82 

Moon  <1925 

Moon  >1925 

6  876 

.52525 

.68170 

.49191 

.68170 

.93 

.70 

.93 

.79 

Many  significant 
for  duplication 

figures  in  weights  are 
of  final  results. 

meaningless 

but  may 

be  required 

Unrealistic.  See  text. 
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TABLE  XII 

Heliocentric  Equatorial  Coordinates.  Mean  Equator  and  Equinox  of  1950.0 


TABLE  XII  (Continued) 
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TABLE  XIII 


Lunar  Osculating  Gee  centric  Equatorial  Elements 


Mean  Equator  and  Equinox  of  1950.0 

Element 

Value  and  Standard 

Deviation 

at  JD  242  0000.5 

a 

.002  563  725  2 

±  .000  000 

e 

.058  227  02 

±  .000  000 

I 

28^680  189 

±°.000  005 

1 

198.587  90 

±  .000  25 

g 

177.267  78 

±  .000  06 

h 

359.000  37 

±  .000  01 

A-8 
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TABLE  XIV 


Osculating  Ecliptic  Elements.  Mean  Ecliptic  and  Equinox  of  1950.0. 


a  [a 

■u.] 

'  e 

I 

1 

g 

It 

Moor 

.002 

563 

725 

2 

.058 

227 

02 

5?2S2 

621 

19^587 

90 

18f.6i6 

29 

354^752 

57 

+.000 

000 

000 

2 

±.000 

000 

09 

+  .000 

005 

±.000 

25 

±.000 

08 

+.000 

06 

.002 

561 

126 

2 

.054 

442 

37 

5.200 

962 

20.647 

11 

166.091 

05 

268.691 

91 

Mercury 

.387 

097 

842 

2 

.205 

624 

5 

7.006 

20 

324.152 

8 

28.836 

5 

47.783 

8 

+.000 

000 

000 

4 

±.000 

000 

3 

±.000 

02 

±.000 

1 

±.000 

2 

±.000 

2 

.387 

098 

432 

3 

.205 

620 

7 

7.002 

53 

355.505 

5 

29.009 

8 

47.707 

7 

Venus 

.723 

325 

568 

1 

.006 

855 

96 

3.394 

451 

271.975 

55 

54.568 

24 

76.331 

82 

±.000 

000 

000 

8 

±.000 

000 

08 

±.000 

007 

±.000 

65 1 

±.000 

66 1 

±.000 

11 

.723 

326 

653 

7 

.006 

796 

79 

3.393 

910 

238.613 

45 

54.968 

21 

76.168 

18 

Earth 

.999 

416 

569 

5 

.016 

945 

64 

.003 

033 

225.239 

00 

92.912 

02 

11.037 

81 

+.000 

000 

000 

8 

±.000 

000 

03 

±.000 

003 

±.000 

11 

±.055 

80* 

±.055 

80' 

.999 

330 

544 

3 

.016 

095 

09 

.002 

544 

309.290 

78 

333.814 

72 

129.478 

95 

Mars 

1.523 

662 

703 

.093 

219 

38 

1.852 

889 

49.430 

12 

285.680 

77 

49.279 

46 

±.000 

000 

002 

±.000 

000 

04 

±.000 

004 

±.000 

02 

±.000 

12 

±.000 

12 

1.523 

634 

008 

.093 

331 

90 

1.848 

036 

57.912 

07 

286.102 

63 

49.098 

38 

Jupiter 

5.202 

965 

95 

.048 

091 

38 

1.307 

499' 

279.423 

14 

273.532 

31 

99.861 

20 

±.000 

000 

03 

±.000 

000 

04 

±.000 

006 

±.000 

05 

±.000 

25 

+.000 

25 

5.202 

854 

11 

.048 

083 

11 

1.306 

451 

307.147 

26 

273.579 

39 

99.981 

55 

Saturn 

9.523 

632 

9 

.053 

680 

17 

2.489 

648 

340.873 

77 

339.048 

64 

1 13.346 

71 

±.000 

000 

2 

±.000 

000 

06 

+.000 

007 

±.000 

06 

±.000 

16 

±.000 

14 

9.528 

437 

1 

.053 

852 

72 

2.490 

769 

358.108 

15 

338.364 

29 

113.178 

83 

Uranus 

19.280 

385 

.044 

255 

51 

.773 

712 

128.117 

5 

100.550 

9 

73.803 

8 

±.000 

003 

±.000 

000 

09 

+.000 

006 

±.000 

2 

±.000 

5 

+.000 

5 

19.172 

780 

.045 

956 

81 

.772 

571 

29.705 

3 

97.248 

0 

73.759 

3 

Neptune  29.985 

79 

.008 

229 

1.775 

952 

89.577 

3 

254.687 

1 

131.229 

2 

+.000 

05 

±.000 

001 

±.000 

006 

±.004 

6  1 

±.004 

8  1 

+.000 

2 

30.129 

69 

.007 

869 

1.773 

438 

220.834 

3 

255.177 

9 

131.300 

2 
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TABLE  XIV  (Continued) 


a  [a.u.] 

e 

i 

1 

8 

h 

Pluto  39.383  9 

.249  779 

177182  67 

248?747 

1 147399 

1097584  06 

±.000  7 

+.000  003 

±.000  04 

±.004 

±.003 

±.000  05 

39.399  3 

.248  065 

17.149  76 

336.543 

114.186 

109.706  18 

Line  1 :  Elements  for  JD  242  0000.5 

Line  2:  Formal  standard  deviations  for  JD  242  0000.5 

Line  3:  Elements  for  JD  244  2000.5 

1)  Small  eccentricity.  Elements  1  and  g  poorly  defined. 

2)  Small  inclination.  Elements  g  and  h  poorly  defined. 
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TABLE  XV 


Correlation  Coefficients 


a 

e 

1 

1 

o 

a 

h 

a 

-.04 

.01 

.03 

.02 

-.03 

c 

-.01 

-.01 

-.03 

.04 

-.01 

Mercury 

I 

.01 

.00 

Svs\ 

331 

-.01 

-.01 

Venus 

1 

.70 

.07 

.02 

^  S 

^  -  .98 

.01 

g 

.00 

-.06 

-.06 

00 

OJ 

1 

-.18 

h 

.00 

.01 

.07 

.00 

-.90  ^ 

a 

.05 

.02 

.39 

.03 

-.02 

e 

.02 

-.01 

-.01 

.01 

.00 

Earth 

I 

.06 

.01 

.01 

.01 

-.01 

Mars 

1 

.17 

.04 

.01 

-.16 

-.01 

g 

.01 

.00 

-.02 

.00 

-.98 

h 

-.01 

.00 

.02 

-.01 

-l.oo" 

a 

-.27 

.00 

-.20 

.13 

.00 

e 

.n" 

.00 

.01 

-.02 

.00 

Jupiter 

I 

.00 

.00 

O 

o 

-.03 

.03 

Saturn 

1 

.07 

.03 

.00 

-.38 

.00 

g 

.01 

.00 

1 

O 

-.21 

-.92 

h 

.01 

.00 

.03 

.00 

-.98 

a 

-.99 

-.01 

-.96 

.96 

-.01 

e 

-.67" 

.01 

.93 

-.93 

.01 

Uranus 

I 

.04 

-.09 

.01 

I 

o 

OO 

.46 

Neptune 

1 

.98 

-.62 

.04 

-1.00 

.01 

g 

-.45 

.25 

.03 

CO 

Tf 

1 

-.05 

h 

.08 

.00 

-.05 

.09 

-.92 

'VSs’nv 

a 

-.85 

-.03 

-.98 

-.55 

-.02 

e 

.79 

..00 

.85 

.39 

.09 

Pluto 

I 

.02 

.02 

^.03 

.44 

-.11 

Moon 

1 

1.00 

.80 

.02 

.37 

.04 

(geocentric) 

g 

-1.00 

-.76 

-.03 

-1.00 

-.24 

equatorial) 

h 

-.04 

.03 

.69 

-  .04 

.02 

Correlations  between  the  six  elements  of  any  one  planet. 
Do  not  read  across  planets. 


TABLE  XVI 


Correlation  Coefficients 


Mercury  Venus  Earth  Mars  Jupiter 


a 

1 

a 

1 

a 

1 

a 

1 

a 

Venus 

a 

.51 

.48 

1.00 

.82 

.14 

.63 

.31 

.13 

Earth 

a 

.58 

.55 

.82 

1.00 

.77 

.37 

.17 

e 

-.10 

.28 

-.23 

1 

.10 

.29 

.14 

.53 

1.00 

.12 

Mars 

a 

.43 

.41 

.63 

.77 

.12 

1.00 

e 

.11 

.13 

.40 

1.00 

1 

.24 

.20 

.31 

-.23 

.37 

-.23 

Jupiter 

C 

a 

-.64 

-.57 

.13 

-.73 

.17 

-.86 

-.14 

-.65 

-.33 

1.00 

-.13 

Moon 

a 

-.68 

-.62 

-.79 

-.93 

-.16 

-.70 

-.36 

-.14 

e 

.58 

.53 

.67 

.79 

.13 

.60 

.31 

.12 

1 

.67 

.61 

.77 

.91 

.15 

.69 

.36 

.14 

2 

.36 

.33 

.42 

.50 

.37 

.20 

-.14 

€1 

© 

-.67 

-.61 

-.78 

-.91 

-.15 

-.69 

-.36 

*2 

-.68 

-.68 

-.62 

-.62 

-.79 

-.79 

-.92 

-.92 

-.15 

-.16 

-.70 

-.70 

-.36 

-.36 

-.14 

-.14 

-.68 

-.62 

-.79 

-.92 

-.15 

-.70 

-.36 

-.14 

6c 

-.67 

-.62 

-.79 

-.92 

-.16 

-.69 

-.36 

-.14 

6c 

-.67 

-.61 

-.78 

-.91 

-.15 

-.69 

-.36 

-.13 

6*7 

-.65 

-.60 

-.76 

-.89 

-.15 

-.67 

-.36 

-.13 

7 

6o 

-.63 

-.58 

-.74 

-.86 

-.15 

-.65 

-.34 

-.13 

e9 

6i  n 

-.58 

-.50 

-.55 

-.47 

-.69 

-.60 

-.80 

-.60 

-.14 

-.12 

-.61 

-.52 

-.32 

-.28 

-.12 

-.10 

M  0 

ell 

-.37 

-.37 

-.45 

-.52 

-.09 

-.40 

-.22 

-.08 
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TABLE  XVII 


Correlation  Coefficients 


Mercury 

Venus 

e 

e 

I 

g 

h 

e 

Venus 

1 

.17 

g 

-.16 

1.00 

Earth 

e 

.31 

.52 

-.28 

1.00 

I 

.24 

-.41 

1 

.17 

-.25 

.52 

g 

-.45 

-.21 

h 

.45 

.21 

Mars 

e 

-.27 

-.13 

-.30 

I 

.23 

1 

-.15 

-.10 

.24 

-.36 

h  .22 


Earth 

I  g 

1.00 

1.00 

.36  -.34 

.35  .34 

-.35  -.34 
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Correlation  Coefficients 


TABLE  XlX 


Atomic  Time 

1912.5  to 

1954.5 

t 

ATf,  — 

Our  Solution 

E.T.0.-U.T.2 

AT-UT 

a 

AT-ET 

1912.5 

1310 

1818 

119 

518 

1913.5 

14.2 

20.1 

5.9 

1914.5 

15.3 

21.4 

6.1 

1915.5 

16.4 

22.7 

6.3 

1916.5 

17.4 

24.0 

1.6 

6.6 

1917.5 

18.3 

24.6 

6.3 

1918.5 

19.1 

25.1 

6.0 

1919.5 

19.8 

25.5 

5.7 

1920.5 

20.5 

25.8 

1.4 

5.3 

1921.5 

21.1 

26.5 

5.4 

1922.5 

21.6 

27.0 

5.4 

1923.5 

22.0 

27.5 

5.5 

1924.5 

22.3 

27.9 

1.1 

5.6 

1925.5 

22.6 

28.1 

5.5 

1926.5 

22.7 

28.1 

5.4 

1927.5 

22.8 

28.0 

5.2 

1928.5 

22.9 

27.9' 

.9 

5.0 

1929.5 

23.0 

27.6 

4.6 

1930.5 

23.2 

27.3 

4.1 

1931.5 

23.3 

27.0 

3.7 

1932.5 

23.5 

26.7 

.7 

3.2 

1933.5 

23.6 

26.5 

2.9 

1934.5 

23.6 

26.3 

2.7 

1935.5 

23.6 

26.1 

2.5 

1936.5 

23.6 

26.0 

.6 

2.4 

1937.5 

23.6 

26.0 

2.4 

1938.5 

23.8 

26.2 

2.4 

1939.5 

24.0 

26.3 

2.3 

1940.5 

24.3 

26.6 

.4 

2.3 

1941.5 

24.7 

26.9 

2.2 

1942.5 

25.2 

27.3 

2.1 

1943.5 

25.6 

27.7 

2.1 

1944.5 

26.1 

28.2 

.3 

2.1 

1945.5 

26.6 

28.6 

2.0 

A-15 
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TABLE  XIX  (Continued) 


t 

AT0  = 

Our  Solution 

E.T.0.-U.T.2 

AT-UT 

a 

AT-ET 

1946.5 

27.1 

29.0 

1.9 

1947.5 

27.6 

29.4 

1.8 

1948.5 

28.2 

29.8 

2 

1.6 

1949.5 

28.9 

30.1 

1.2 

1950.5 

29.4 

30.3 

.9 

1951.5 

29.7 

30.6 

.9 

1952.5 

30.3 

30.8 

.2 

.5 

1953.5 

31.0 

30.9 

-.1 

1954.5 

31.1 

31.1 

.0 
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APPENDIX  B 
AUXILIARY  FIGURES 
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Moon.  Residuals  in  Right  Ascension  Times  cos 


Residuals  in  Right  Ascension  Times  cos 


B-4 
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/ 


Mercury.  Residuals  in  Right  Ascension 


Venus.  Residuals  in  Right  Ascension  Times  cos  5 


Residuals  in  Right  Ascension  Times  cos 


Mars.  Residuals  in  Declination 
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Jupiter.  Residuals  in  Declination 


Saturn.  Residuals  in  Right  Ascension  Times  cos 


Uranus.  Residuals  in  Right  Ascension  Times  cos 


Neptune.  Residuals  in  Right  Ascension  Times  cos 
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Pluto.  Residuals  in  Right  Ascension  Times  cos 
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Pluto.  Residuals  in  Declination 


